-
Notifications
You must be signed in to change notification settings - Fork 0
/
report.Rnw
1529 lines (1209 loc) · 68.9 KB
/
report.Rnw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
\documentclass{article}
% skomentowana dla MAC wersja {inputenc}
\usepackage[utf8]{inputenc}
% \usepackage[T1]{fontenc}
%\usepackage[cp1250]{inputenc}
\usepackage{amssymb}
\usepackage{color}
\usepackage{amsmath}
\usepackage{Sweave}
\usepackage{enumerate}
\usepackage{hyperref}
\usepackage{polski}
\title{California Housing: Statystyczna analiza danych}
\author{\textbf{404811, Mykola Haltiuk}, czwartek $17^{50}$\\
\textit{AGH, Wydział Informatyki Elektroniki i Telekomunikacji}\\
\textit{Rachunek prawdopodobieństwa i statystyka 2020/2021}}
\date{Kraków, \today}
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
\textit{Ja, niżej podpisany(na) własnoręcznym podpisem deklaruję, że przygotowałem(łam) przedstawiony do oceny projekt samodzielnie i żadna jego część nie jest kopią pracy innej osoby.}
\begin{flushright}
{............................................}
\end{flushright}
\section{Streszczenie raportu}
Raport powstał w oparciu o analizę danych dotyczących cen na nieruchomoś\'c w
Kalifornii stanem na 1990 rok.
\section{Opis danych}
Dane do projektu pochodzą ze strony \href{https://www.kaggle.com/camnugent/california-housing-prices}{\texttt{Kaggle}}. Są one przedstawione w pliku \textit{.csv}. Ze strony można
spokojnie ich pobra\'c.
\noindent
\quad Otóż, dane ska\l adają si\c ez 20640 rekordów o 10 cechach. Każdy rekord
odpowiada za pewny blok domów. Opis poszczególnych cech:
\begin{itemize}
\item \textbf{longitude} - d\l ugoś\'c geograficzna bloku (\textit{numeryczna})
\item \textbf{latitude} - szerokoś\'c geograficzna bloku (\textit{numeryczna})
\item \textbf{housing\_median\_age} - mediana wieków domów w bloku (\textit{numeryczna})
\item \textbf{total\_rooms} - iloś\'c pokoi w bloku (\textit{numeryczna})
\item \textbf{total\_bedrooms} - iloś\'c sypaleń w bloku (\textit{numeryczna})
\item \textbf{population} - populacja w bloku (iloś\'c osób) (\textit{numeryczna})
\item \textbf{households} - iloś\'c gospodarstw w bloku (\textit{numeryczna})
\item \textbf{median\_income} - mediana zarobków w tym bloku (\textit{numeryczna})
\item \textbf{median\_house\_value} - mediana cen na dom (\textit{numeryczna})
\item \textbf{ocean\_proximity} - blizoś\'c do oceanu (\textit{kategoryczna})
\end{itemize}
\noindent
\quad Dane nie są gotowe do analizy, potrzebują dużo modyfikowania i czyszczenia. Bardzo dużo jest wartości odstających, praktycznie wszystkie są z góry, czyli wi\c eksze od mediany, co psuje rozk\l ady, które już nie da si\c e nazwa\'c normalnymi (b\c edzie omówiono później i bardziej szczegó\l owo).
\noindent
\quad Istnieje kilka problemów z cechami:
\begin{itemize}
\item \textbf{housing\_median\_age} jak i \textbf{median\_house\_value} są ograniczone od góry (\textit{capped}). Czyli wszystkie wartości powyżej pewnej wartości $x$ są zmniejszone do tej wartości. To powoduje duże skoncentrowanie w pewnych punktach.
\item \textbf{total\_bedrooms} ma niektóre wartości NA, czyli są nieznane (\textit{missing values}).
\item \textbf{ocean\_proximity} jest zmienną kategoryczną o 5 różnych wartości, które nie są uporządkowane, co może sprawi\'c problem, bo trudno b\c edzie przekszta\l ci\'c tą zmienną kategoryczną na numeryczną.
\end{itemize}
\noindent
\quad Samo przygotowanie danych zostanie opisane troszeczk\c e niżej w sekcji \textbf{Data cleaning}.
\noindent
\quad Także niektóre cechy nie pokazują ciekawej informacji, wi\c ec zostaną reogranizowane i usuni\c ete, a na ich miejscu zostaną stworzone nowe. Czyli zrobimy tak zwany \textit{feature engineering}. Sam proces zostanie też opisany niżej.
\noindent
\quad Przed przystąpieniem do przygotowywania danych spróbujmy ich za\l adowa\'c i zobaczy\'c jak wyglądają.
<<echo=FALSE>>=
library(corrplot)
library(purrr)
library(ggplot2)
library(hrbrthemes)
library(raster)
library(tidyverse)
library(viridis)
@
<<loading_data>>=
fpath <- "housing.csv"
all_housing <- read.csv(fpath,header=TRUE,stringsAsFactors = FALSE)
@
\noindent
\quad Otóż, plik przez nas za\l adowany ma nazw\c e \textit{housing.csv}, wczytujemy go jako \textbf{data.frame} zachowując nazwy kolumn.
\noindent
\quad Zobaczymy czy zgadzają si\c e wymiary z przedstawionymi liczbami wyżej.
<<dim_all>>=
dim(all_housing)
@
\noindent
\quad Jak widzimy, to rzeczywiście: mamy 20640 rekordów o 10 cechach, czyli 20640 wierszy i 10 kolumn. Zobaczymy także czy nazwy odpowiadają przedstawionym wyżej:
<<colnames_all>>=
colnames(all_housing)
@
\noindent
\quad Wszystko zgadza si\c e. W końcu zobaczymy ma\l y opis tych danych i 4 pierwsze rekordy, żeby dobrze rozumie\'c czym są i jak wyglądają.
<<str_all>>=
str(all_housing)
@
<<head_all>>=
head(all_housing, 4)
@
\noindent
\quad Wszystkie kolumny, oprócz \textbf{ocean\_proximity} są numeryczne, ostatnia jest cechą kategoryczną, wi\c ec musimy to później naprawi\'c. W tej chwili możemy sobie zobaczy\'c wszystkie możliwe wartości tej zmiennej:
<<levels_all>>=
levels(as.factor(all_housing$ocean_proximity))
@
\noindent
\quad Otóż, mamy 5 możliwych wartości:
\begin{itemize}
\item \textbf{INLAND} - w środku Kalifornii, daleko od oceanu
\item \textbf{<1H OCEAN} - w jednej godzinie jazdy od oceanu
\item \textbf{NEAR OCEAN} - blisko oceanu
\item \textbf{NEAR BAY} - blisko zatoki (San-Francisco)
\item \textbf{ISLAND} - wyspa
\end{itemize}
\noindent
\quad Dok\l adniej i na karcie to zostanie przedstawione poźniej, gdy b\c edziemy zajmowa\'c si\c e tą zmienną kategoryczną.
\noindent
\quad Otóż, otrzymaliśmy jakąś pierwotną informacj\c e na temat tego zbioru. Niżej spróbujemy go przeanalizowa\'c jeszcze g\l\c ebiej, żeby przygotowa\'c do prawdziwej statystycznej analizy.
\section{Data cleaning}
\subsection{\textit{NA} wartości}
\quad Sprawdźmy sobie ile jest nieznanych wartości w każdej kolumnie.
<<na_columns_all>>=
colSums(is.na(all_housing))
@
\noindent
\quad Jak i by\l o napisane w liście problemów tego zbioru, w kolumnie \textbf{total\_bedrooms} pojawia si\c e 207 nieznanych wartości, z którymi musimy coś zrobi\' c. Istnieją różne sposoby na rozwiązanie tego problemu:
\begin{itemize}
\item Usuni\c ecie wszystkich rekordów posiadających nieznane wartości. W naszym przypadku to jest ca\l kiem sensowne podejście, bo takich rekordów jest dużo mniej w stosunku do ilości wszystkich.
\item Usuni\c ecie ca\l ej kolumny, ale tak tracimy dużo ciekawej informacji, wi\c e to odrzucamy od razu.
\item Przypisanie jakiejś wartości w miejscu nieznanych wartości, np. 0. To nie jest dobrym podejściem i bardzo zależy od samych danych, bo jeżeli dane są z przedzia\l u $[1200, 1250]$, to otrzymujemy bardzo dużo wartości odstających, które psują nasz rozk\l ad. Tą możliwoś\'c od razu odrzucamy.
\item Przypisanie jakiejś statystyki rozk\l adu w miejscu nieznanych wartości, czyli mediany, wartości średniej i tp. To jest bardzo sensowne podejście i stosuje si\c e najcz\c eściej.
\item Analiza zależności z innymi cechami i sprawdzenie możliwości przewidywania nieznanej wartości na podstawie innych cech. To podejście jest rzadko stosowane, ale jest ciekawe i pozwala nam doś\'c logicznie podebra\'c te wartości, a nie przypisa\'c po prostu coś.
\end{itemize}
\noindent
\quad Otóż, o ile chcia\l oby si\c e zachowa\'c wszystkie rekordy, to odrzucamy pierwszą możliwoś\'c. O ile chcemy przeprowadza\'c testy, to 200 jednakowych wartości może troch\c e wp\l ywa\'c na wynik, i o ile dodatkowo ostatnie podejście jest najciekawsze, to spróbujemy stworzy\'c minimalny model regresyjny, żeby podebra\'c nieznane wartości, ale na początku musimy sprawdzi\'c czy zdążymy to zrobi\'c bazując si\c e tylko na jednej zmiennej, czyli chcemy sprawdzi\'c jak silnie wp\l ywają inne cechy na \textbf{total\_bedrooms} i czy są silnie skorelowane.
\noindent
\quad Żeby policzy\'c macierz wspó\l czynników korelacji, musimy na początku zostawi\'c tylko numeryczne cechy, a przed tym jeszcze usuną\'c rekordy z nieznanymi wartościami:
<<non-na_only_numeric, fig=TRUE>>=
cleaned <- all_housing[rowSums(is.na(all_housing)) == 0,]
colSums(is.na(cleaned)) # rechecking
cleanedNumeric <- cleaned[ , purrr::map_lgl(cleaned, is.numeric)]
# calculating the correlation matrix
hcor <- cor(cleanedNumeric)
# plotting
corrplot(hcor, method='number')
@
\noindent
\quad Najbardziej nas interesuje kolumna/wiersz \textbf{total\_bedrooms}, dlatego, że musimy znaleź\'c inną cech\c e, która jest bardzo skorelowana z \textbf{total\_bedrooms}. Jak widzimy, cecha \textbf{households} ma wspó\l czynnik korelacji z interesującą nas kolumną równy 0.98, co jest bardzo dużą wartością i możemy z pewnością powiedzie\'c, że z doś\'c silnym prawdopodobieństwiem zgadniemy wartoś\'c bliską rzeczywistej \textbf{total\_bedrooms} mając wartoś\'c \textbf{households}. Wspó\l czynnik korelacji bliski 1 wskazuje, że im wi\c ecej gospodarstw, tym wi\c ecej sypaleń, co wydaje si\c e by\'c doś\'c logiczne.
\noindent
\quad Warto także zauważy\'c, że mamy dużo cech silnie skorelowanych, co nie wp\l ywa dobrze na tworzenie różnych modeli regresyjnych i klasyfikacyjnych, bo zwi\c eksza wymiarowoś\'c rekordów nie zwi\c ekszając efektywnej informacji. Ten problem rozwiążemy w \textbf{Feature engineering}.
\noindent
\quad Wró\'cmy do naszych nieznanych wartości. Spróbujmy narysowa\'c sobie wykres zależności tych dwóch cech i dodatkowo sprawdźmy dok\l adniejszą wartoś\'c wspó\l czynnika korelacji.
<<brcor>>=
brcor <- cor(cleaned$total_bedrooms, cleaned$households)
brcor
@
\begin{figure}[h!]
\centering
<<bed_house_dependence, fig=TRUE, include=FALSE, width=8, height=6>>=
ggplot(cleaned, aes(x=households, y=total_bedrooms)) +
labs(title="Total Bedrooms / Households Dependence",
x="Households", y="Total Bedrooms") +
theme(plot.title=element_text(size=30, face="bold", family='sans'),
axis.text.x=element_text(size=15, family='sans'),
axis.text.y=element_text(size=15, family='sans'),
axis.title.x=element_text(size=25, family='sans'),
axis.title.y=element_text(size=25, family='sans')) +
geom_point( color="#69b3a2" ) +
theme_ipsum(base_family = 'sans')
@
\includegraphics[width=\textwidth]{report-bed_house_dependence}
\end{figure}
\noindent
\quad Otóż, widzimy, że zależnoś\'c jest naprawd\c e bliska liniowej. \textit{Ma\l y komentarz: taki d\l ugi kod zosta\l\ wykorzystany do rysowania wi\c ekszości wykresów, wi\c ec dalej nie b\c edzie przedstawiony, a b\c edą pokazane same wykresy, żeby nie zaśmieca\'c powtarzającym si\c e kodem ten raport}.
\noindent
\quad Zgodnie z zapowiedzią, tworzymy prosty model regresyjny:
<<regressNA>>=
bhmodel <- lm(total_bedrooms ~ households, cleaned)
bhmodel
@
\noindent
\quad Narysujmy sobie to od razu z prostą stworzoną podczas regresji:
\begin{figure}[h!]
\centering
<<regfig, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
ggplot(cleaned, aes(x=households, y=total_bedrooms)) +
labs(title="Total Bedrooms / Households Dependence",
x="Households", y="Total Bedrooms") +
theme(plot.title=element_text(size=30, face="bold", family='sans'),
axis.text.x=element_text(size=15, family='sans'),
axis.text.y=element_text(size=15, family='sans'),
axis.title.x=element_text(size=25, family='sans'),
axis.title.y=element_text(size=25, family='sans')) +
geom_point( color="#69b3a2" ) +
geom_smooth(method=lm , color="darkorange", se=FALSE) +
theme_ipsum(base_family = 'sans')
@
\includegraphics[width=\textwidth]{report-regfig}
\end{figure}
\noindent
\quad Model zosta\l\ stworzony, wi\c ec zosta\l o si\c e tylko zamieni\'c wszystkie \textit{NA} na wartości otrzymane z modelu.
<<substNA>>=
# replacing the missing values with predictions of our model
all_housing$total_bedrooms[is.na(all_housing$total_bedrooms)] =
predict(bhmodel, all_housing[is.na(all_housing$total_bedrooms), ])
# checking if the null values has left
colSums(is.na(all_housing))
@
\subsection{Cecha kategoryczna \textit{ocean\_proximity}}
\quad Pora zają\'c si\c e \textbf{ocean\_proximity}, zmienną kategoryczną, którą musimy przekszta\l ci\'c na numeryczną. Istnieją na to różne sposoby:
\begin{itemize}
\item \textbf{Integer encoding} - zamiana wartości kategorycznych na numeryczne wprost stosując pewną funkcj\c e mapującą. To podejście jest stosowane, gdy jest oczywisty pewny porządek wartości kategorycznych (np. \textit{low} - \textit{medium} - \textit{high}).
\item \textbf{One-hot encoding} - dekompozycja jednej cechy na kilka tak, żeby każda możliwa wartoś\'c kategoryczna przedstawia\l a nową kolumn\c e. I wtedy dla każdej kolumny przypisujemy \textit{bool}-owe wartości, czyli pewien rekord mia\l\ tą wartoś\'c kategoryczną lub nie mia\l.
\end{itemize}
\noindent
\quad Chociaż nasz przypadek bardziej odpowiada podejściu one-hot, nie b\c edziemy komplikowa\'c naszych danych, a po prostu spróbujmy stworzy\'c taki porządek możliwych wartości \textbf{ocean\_proximity}, żeby można by\l o go uzasadni\'c, żeby nie by\l\ przypadkowym.
\noindent
\quad Na początku mając d\l ugoś\'c i szerokoś\'c geograficzną narysujmy sobie jak są rozrzucone te rekordy z różnymi wartościami \textbf{ocean\_proximity}.
\begin{figure}[h!]
\centering
<<califmapproxim, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
states <- c('California')
us <- getData("GADM",country="USA",level=1)
us.states <- us[us$NAME_1 %in% states,]
us.bbox <- bbox(us.states)
xlim <- c(min(us.bbox[1,1]),max(us.bbox[1,2]))
ylim <- c(min(us.bbox[2,1]),max(us.bbox[2,2]))
plot(us.states, xlim=xlim, ylim=ylim, family="sans", main="Ocean proximity")
points(all_housing$longitude[all_housing$ocean_proximity=="<1H OCEAN"],
all_housing$latitude[all_housing$ocean_proximity=="<1H OCEAN"],
col = "#F8756D", cex = .5)
points(all_housing$longitude[all_housing$ocean_proximity=="NEAR BAY"],
all_housing$latitude[all_housing$ocean_proximity=="NEAR BAY"],
col = "#00BE7D", cex = .5)
points(all_housing$longitude[all_housing$ocean_proximity=="INLAND"],
all_housing$latitude[all_housing$ocean_proximity=="INLAND"],
col = "#A3A500", cex = .5)
points(all_housing$longitude[all_housing$ocean_proximity=="ISLAND"],
all_housing$latitude[all_housing$ocean_proximity=="ISLAND"],
col = "#03B0F6", cex = .5)
points(all_housing$longitude[all_housing$ocean_proximity=="NEAR OCEAN"],
all_housing$latitude[all_housing$ocean_proximity=="NEAR OCEAN"],
col = "#E76AF3", cex = .5)
legend(legend = c("<1H OCEAN", "NEAR BAY", "INLAND", "ISLAND", "NEAR OCEAN"),
fill = c("#F8756D", "#00BE7D", "#A3A500", "#03B0F6", "#E76AF3"),
"topright")
@
\includegraphics[width=0.68\textwidth]{report-califmapproxim}
\end{figure}
\noindent
\quad Tworzymy funkcj\c e mapującą w oparciu o \textbf{median\_house\_value}, bo w\l aśnie cel tego zbioru danych - stworzenie modelu do przewidywania mediany cen domów, żeby firma-klient mog\l a zdecydowa\'c si\c e lub nie na inwestycje.
\noindent
\quad Otóż, początkowo bazując si\c e na w\l asnym odczuciu stworzymy taki porządek: \textit{INLAND} $\rightarrow$ \textit{<1H OCEAN} $\rightarrow$ \textit{NEAR BAY} $\rightarrow$ \textit{NEAR OCEAN} $\rightarrow$ \textit{ISLAND}. Uzasadnienie temu porządkowi jest to, że im bliżej do oceanu, tym droższe stają si\c e mieszkania.
\noindent
\quad Ale nie możemy budowa\'c analizy na w\l asnych odczuciach, wi\c ec musimy sprawdzic\'c jakoś naszą propozycj\c e uporządkowania. Porównajmy mediany cen dla poszczególnych wartości kategorycznych i narysujmy od razu ca\l e \textit{box-plot}'y dla porównania:
\begin{figure}[h!]
\centering
<<proximityboxes, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
all_housing %>% ggplot(aes(x=ocean_proximity, y=median_house_value,
fill=ocean_proximity)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha=0.6) +
geom_jitter(color="darkorange", size=0.4, alpha=0.4) +
theme_ipsum(base_family = 'sans') +
theme(
legend.position="none",
plot.title = element_text(size=20),
axis.text.x=element_text(size=12),
axis.text.y=element_text(size=12),
axis.title.x=element_text(size=15),
axis.title.y=element_text(size=15)
) +
labs(title="Median House Value by the Ocean Proximity",
x="Ocean Proximity", y='Median House Value')
@
\includegraphics[width=\textwidth]{report-proximityboxes}
\end{figure}
\noindent
\quad Byliśmy doś\'c blisko. Okaza\l o si\c e, że \textit{NEAR BAY} ma median\c e wi\c ekszą od mediany \textit{NEAR OCEAN}. Pierwsze i trzecie kwantyle też odpowiadają naszemu porządkowi. Przewaga \textit{NEAR BAY} może by\'c objaśniona urbanizacją brzegu zatoki. Takie miasta jak San Francisco, Berkley, Oakland, Richmond, Fremont, San Jose, Silicon Valley są rozrzucone po brzegu zatoki. Sacramento też jest doś\'c blisko, a to powoduje wzrost cen na mieszkanie.
\noindent
\quad Otóż, nasza funkcja mapująca wygląda nast\c epująco: \textit{INLAND} $\rightarrow$ 1, \textit{<1H OCEAN} $\rightarrow$ 2, \textit{NEAR OCEAN} $\rightarrow$ 3, \textit{NEAR BAY} $\rightarrow$ 4, \textit{ISLAND} $\rightarrow$ 5.
\noindent
\quad Zamieniamy naszą cech\c e kategoryczną na numeryczną.
<<proxtonum>>=
fac <- factor(all_housing$ocean_proximity,
levels=c('INLAND', '<1H OCEAN', 'NEAR OCEAN', 'NEAR BAY', 'ISLAND'))
# let's check the mapping
data.frame(levels = unique(fac), value = as.numeric(unique(fac)))
all_housing$ocean_proximity <- as.numeric(fac)
@
\subsection{Wartości odstające (\textit{outliers})}
\quad Ze wszystkich wymienionych na początku problemów zosta\l\ si\c e tylko jeden - ograniczone wartości (\textit{capped}). Z box-plotu w poprzedniej sekcji możemy zobaczy\'c jak dużo wartości cen są na poziomie 500 tysi\c ecy dolarów. Jest to dlatego, jak by\l o wcześniej już powiedziane, że wszystkie ceny powyżej tego progu zosta\l y zapisane równymi temu progowi. To tworzy pewne problemy z rozk\l adem i może wp\l yną\'c na regresj\c e doś\'c znacząco.
\noindent
\quad Znajdźmy ten próg i sprawdźmy ile wartości mu odpowiadają:
<<cappedvals>>=
lim <- max(all_housing$median_house_value)
lim
sum(all_housing$median_house_value == lim)
@
\noindent
\quad Czyli tak naprawd\c e, jeżeli rzeczywisty rozk\l ad wartości jest $X$, to nasz $Y = min(X, 500k)$.
\noindent
\quad Co możemy z tym zrobi\'c?
\begin{itemize}
\item Możemy ich usuną\'c tracąc 5\% danych i zamiast dużej ściany rekordów w tym miejscu powstanie przepaś\'c, co też nie jest najlepszym rozwiązaniem.
\item W przypadku gdy nie potrzebujemy, żeby model przewidywa\l\ wartości powyżej 500k, to najlepszym rozwiązaniem b\c edzie zostawienie tego jak jest.
\item W innym przypadku, musielibyśmy znaleź\'c prawdziwe wartości dla tych rekordów co dla nas nie jest ani istotne, ani możliwe, wi\c ec tą opcj\c e od razu odrzucamy.
\end{itemize}
\noindent
\quad Usuwamy czy zostawiamy? Dla potrzeb statystycznych zostawimy sobie te dane, bo mogą by\'c ciekawe później, zawsze b\c edziemy mieli możliwoś\'c odfiltrowania tych outlier'ów.
\noindent
\quad Także zauważamy, że pomys\l\ porównywa\'c wartości kategoryczne przy pomocy mediany by\l\ dobrym, bo wartoś\'c średnia psuje si\c e przy takiej zamianie rozk\l adów i ograniczaniu wartości, a mediana zostaje si\c e taką samą dopóki ta zmiana rozk\l adu nie zmienia po\l ow\c e wartości.
\noindent
\quad Otóż, wszystkie problemy zosta\l y rozwiązane lub przynajmniej rozpatrzone, nasz zbiór danych już jest zmodyfikowany. Możemy go sobie zapisa\'c jako odr\c ebny plik, żeby potem nie powtarza\'c tego samego.
<<writingcleaned, eval=FALSE>>=
write.csv(all_housing, file="cleaned_housing.csv", sep=",",
col.names=TRUE, row.names=FALSE)
@
\section{Feature engineering}
\subsection{Tworzenie nowych cech}
<<echo=FALSE>>=
library(corrplot)
library(purrr)
library(ggplot2)
library(hrbrthemes)
library(raster)
library(tidyverse)
library(viridis)
library(forcats)
library(dplyr)
library(ggpubr)
housing <- read.csv("cleaned_housing.csv", header=TRUE)
density.plot.properties <- list(
geom_density(fill="darkorange", color="darkorange", alpha=1),
theme_ipsum_rc(grid='Y', base_family = 'sans'),
theme(
legend.position="none",
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(size=5),
axis.text.y = element_text(size=8)
))
@
\noindent
\quad Otóż, spróbujmy sobie narysowa\'c rozk\l ady wszystkich cech:
\begin{figure}[h!]
\centering
<<distrs, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 8>>=
# longitude
pb1 <- housing %>%
ggplot( aes(x=longitude)) + density.plot.properties
# latitude
pb2 <- housing %>%
ggplot( aes(x=latitude)) + density.plot.properties
# housing_median_age
pb3 <- housing %>%
ggplot( aes(x=housing_median_age)) + density.plot.properties
# total_rooms
pb4 <- housing %>%
ggplot( aes(x=total_rooms)) + density.plot.properties
pb5 <- housing %>%
filter(total_rooms < 10000) %>%
ggplot( aes(x=total_rooms)) + density.plot.properties
# total_bedrooms
pb6 <- housing %>%
ggplot( aes(x=total_bedrooms)) + density.plot.properties
pb7 <- housing %>%
filter(total_bedrooms < 2000) %>%
ggplot( aes(x=total_bedrooms)) + density.plot.properties
# population
pb8 <- housing %>%
ggplot( aes(x=population)) + density.plot.properties
pb9 <- housing %>%
filter(population < 5000) %>%
ggplot( aes(x=population)) + density.plot.properties
# households
pb10 <- housing %>%
ggplot( aes(x=households)) + density.plot.properties
pb11 <- housing %>%
filter(households < 2000) %>%
ggplot( aes(x=households)) + density.plot.properties
# median_income
pb12 <- housing %>%
ggplot( aes(x=median_income)) + density.plot.properties
# median_house_value
pb13 <- housing %>%
ggplot( aes(x=median_house_value)) + density.plot.properties
# ocean_proximity
pb14 <- housing %>%
ggplot( aes(x=ocean_proximity)) + density.plot.properties
# drawing it all at one page
theme_set(theme_pubr(base_family = 'sans'))
ggarrange(pb1, pb2, pb3, pb4, pb5, pb6, pb7, pb8, pb9, pb10, pb11, pb12, pb13,
pb14, ncol=4, nrow=4, font.label=list(size = 8, family='sans'),
labels=c('Longitude', 'Latitude', 'Median Age', 'Total Rooms',
'Total Rooms filtered (< 10k)', 'Total Bedrooms',
'Total Bedrooms filtered (< 2k)', 'Population',
'Population filtered (< 5k)', 'Households',
'Households filtered (< 2k)', 'Median Income',
'Median House Value', 'Ocean Proximity'))
@
\includegraphics[width=0.98\textwidth]{report-distrs}
\end{figure}
\noindent
\quad Jak widzimy, dużo rozk\l adów mają ci\c eżki prawy ogon, czyli są \textit{right-tailed}, \textit{heavy-tailed}. Musimy coś z tym zrobi\'c jeżeli chcemy efektywnie testowa\'c i stosowa\'c dane z tych kolumn/cech. Także wcześniej zauważyliśmy, że \textbf{median\_house\_value} jest bardzo s\l abo powiązana z wi\c eksząścią cech, co nie u\l atwia nam życia. Wi\c ec spróbujmy rozwiąza\'c od razu kilka problemów przy pomocy \textit{feature engineering}, czyli tworzenie/modyfikacja/usuwanie kolumn tak, żeby dojś\'c do stanu z lepszą wygodą dla nas.
\noindent
\quad Patrząc na podane nam kolumny, możemy śmia\l o stwierdzi\'c, że \textbf{total\_rooms}, \textbf{total\_bedrooms} absolutnie nic nie wnoszą z praktycznego punktu widzenia, ponieważ bardzo zależą od \textbf{households}. Te duże wartości wspó\l czynników korelacji też naprawimy w ten sposób. Ótóż, przetestujmy sobie takie nowe kolumny/cechy:
\begin{itemize}
\item \textbf{bedrooms\_per\_rooms} - iloś\'c sypalni na pokój
\item \textbf{rooms\_per\_household} - iloś\'c pokoji na gospodarstwo
\item \textbf{bedrooms\_per\_person} - iloś\'c sypalni na osob\c e
\item \textbf{rooms\_per\_person} - iloś\'c pokoji na osob\c e
\item \textbf{bedrooms\_per\_household} - iloś\'c sypalni na gospodarstwo
\item \textbf{people\_per\_household} - iloś\'c osób na gospodarstwo
\end{itemize}
\noindent
\quad Oczywiście, że nie możemy zostawi\'c wszystkie cechy, musimy usuną\'c ich tyle, żeby nie straci\'c informacji i jednocześnie jakaś cecha nie by\l a zdefiniowana jako liniowa kombinacja innych cech.
\noindent
\quad Moglibyśmy stworzy\'c jeszcze i inne cechy. Dobrym pomys\l em by\l oby stworzenie cechy wyznaczającej odleg\l oś\'c do najbliższego dużego miasta i podobne geograficzne cechy, bo możemy to zrobi\'c stosując \textbf{longitude} i \textbf{latitude}. Ale to jest już nadmiar dla projektu naszej skali, wi\c ec kontynuujmy z danymi, które mamy.
\noindent
\quad Otóż, stwórzmy takie nowe dane:
<<newf>>=
housing$bedrooms_per_rooms <- housing$total_bedrooms / housing$total_rooms
housing$rooms_per_household <- housing$total_rooms / housing$households
housing$bedrooms_per_person <- housing$total_bedrooms / housing$population
housing$rooms_per_person <- housing$total_rooms / housing$population
housing$bedrooms_per_household <- housing$total_bedrooms / housing$households
housing$people_per_household <- housing$population / housing$households
@
\newpage
\noindent
\quad Sprawdźmy teraz macierz wspó\l czynników korelacji pomi\c edzy nowymi i starymi cechami:
\begin{figure}[h!]
\centering
<<newfcorr, echo=FALSE, fig=TRUE, include=FALSE, width = 12, height = 9>>=
corfe <- cor(housing)
corrplot(corfe, method='number')
@
\includegraphics[width=0.9\textwidth]{report-newfcorr}
\end{figure}
\noindent
\quad Także sprawdźmy jak wyglądają rozk\l ady nowych cech:
\begin{figure}[h!]
\centering
<<newfdistrs, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
# bedrooms_per_rooms
pfe1 <- housing %>%
ggplot( aes(x=bedrooms_per_rooms)) + density.plot.properties
# rooms_per_household
pfe2 <- housing %>%
ggplot( aes(x=rooms_per_household)) + density.plot.properties
# bedrooms_per_person
pfe3 <- housing %>%
ggplot( aes(x=bedrooms_per_person)) + density.plot.properties
# rooms_per_person
pfe4 <- housing %>%
ggplot( aes(x=rooms_per_person)) + density.plot.properties
# people_per_household
pfe5 <- housing %>%
ggplot( aes(x=people_per_household)) + density.plot.properties
# bedrooms_per_household
pfe6 <- housing %>%
ggplot( aes(x=bedrooms_per_household)) + density.plot.properties
# drawing it all at one page
ggarrange(pfe1, pfe2, pfe3, pfe4, pfe5, pfe6,
ncol=2, nrow=3, font.label=list(size = 12, family='sans'),
labels=c('Bedrooms per room', 'Rooms per household',
'Bedrooms per person', 'Rooms per person',
'People per household', 'Bedrooms per household'))
@
\includegraphics[width=0.9\textwidth]{report-newfdistrs}
\end{figure}
\noindent
\quad Widzimy, że jest dużo outlier'ów, wi\c ec na początku odfiltrowa\'c tak, żeby zobaczy\'c bliżej najwi\c ekszą koncentracj\c e punktów.
\begin{figure}[h!]
\centering
<<newfdistrsfilt, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
# bedrooms_per_rooms
pfe1f <- housing %>%
filter(bedrooms_per_rooms < 0.4) %>%
ggplot( aes(x=bedrooms_per_rooms)) + density.plot.properties
# rooms_per_household
pfe2f <- housing %>%
filter(rooms_per_household < 11) %>%
ggplot( aes(x=rooms_per_household)) + density.plot.properties
# bedrooms_per_person
pfe3f <- housing %>%
filter(bedrooms_per_person < 1) %>%
ggplot( aes(x=bedrooms_per_person)) + density.plot.properties
# rooms_per_person
pfe4f <- housing %>%
filter(rooms_per_person < 4) %>%
ggplot( aes(x=rooms_per_person)) + density.plot.properties
# people_per_household
pfe5f <- housing %>%
filter(people_per_household < 6) %>%
ggplot( aes(x=people_per_household)) + density.plot.properties
# bedrooms_per_household
pfe6f <- housing %>%
filter(bedrooms_per_household < 2) %>%
ggplot( aes(x=bedrooms_per_household)) + density.plot.properties
# drawing it all at one page
ggarrange(pfe1f, pfe2f, pfe3f, pfe4f, pfe5f, pfe6f,
ncol=2, nrow=3, font.label=list(size = 12, family='sans'),
labels=c('Bedrooms per room filtered', 'Rooms per household filtered',
'Bedrooms per person filtered', 'Rooms per person filtered',
'People per household filtered',
'Bedrooms per household filtered'))
@
\includegraphics[width=\textwidth]{report-newfdistrsfilt}
\end{figure}
\noindent
\quad Otóż, musimy sobie wybra\'c cechy, które zostawiamy, a które zostaną wyrzucone.
\noindent
\quad Jak widzimy, \textbf{bedrooms\_per\_person} nic nie wnoszą, wspó\l czynnik korelacji z \textbf{median\_house\_value} jest bardzo niski, wi\c ec wyrzucamy tą cech\c e. Dużo ciekawszym jest dla nas \textbf{bedrooms\_per\_rooms}, bo ma doś\'c znaczący wsp. korelacji z interesującą nas cechą, a także nie ma bliskich do liniowych zależności z innymi cechami. Wi\c ec usuwamy wszystkie cechy powiązane z sypalniami, oprócz tej.
\noindent
\quad Także ciekawą dla nas cechą jest \textbf{rooms\_per\_person}, bo ma wsp. korelacji równy 0.21, i jest bardziej wygodną, niż \textbf{total\_rooms}. Dodatkowo pokazuje ciekawsze dane. Zmieniamy \textbf{total\_rooms} na nią.
\noindent
\quad Inną ciekawą cechą nie w jakości znaczącej dla regresji, tylko tym, że ma niskie wsp. korelacji z cechami, bazując si\c e na których, przewidujemy wartoś\'c \textbf{median\_house\_value}, także tym, że ma doś\'c bliski do normalnego rozk\l ad, przynajmniej tak to wygląda na wykresie.
\noindent
\quad Reszt\c e niepotrzebnych cech usuwamy. Popatrzmy jak wygląda w tej chwili nasz zbiór danych.
<<feended>>=
housingf <- dplyr::select(housing, -c('total_rooms', 'population', 'total_bedrooms',
'bedrooms_per_household', 'rooms_per_household',
'bedrooms_per_person'))
str(housingf)
@
\noindent
\quad Także zobaczmy macierz wsp. korelacji po \textit{feature engineering}.
\begin{figure}[h!]
\centering
<<feendedcorr, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
fedone <- cor(housingf)
corrplot(fedone, method='number')
@
\includegraphics[width=\textwidth]{report-feendedcorr}
\end{figure}
\noindent
\quad Podsumowując, możemy stwierzdi\'c, że pokonaliśmy dobrą robot\c e: nie mamy cech, które są silnie powiązane, oprócz \textbf{longitude} i \textbf{latitude}, ale to wynika z formy Kalifornii, i \textbf{bedrooms\_per\_rooms} i \textbf{median\_income}, ale to nie wp\l ynie doś\'c znacząco na nasz model, wi\c ec możemy ich zostawi\'c. W końcu, mamy dużo nowych cech o dużo lepszej korelacji z \textbf{median\_house\_value}.
\noindent
\quad Zapisujemy zmodyfikowane dane do kolejnego pliku.
<<fewritecsv, eval=FALSE>>=
write.csv(housingf, file="fe_housing.csv", sep=",",
col.names=TRUE, row.names=FALSE)
@
\subsection{Dodatkowe możliwy zmiany}
\quad Pozostają nam kilka problemów, które zauważyliśmy:
\begin{itemize}
\item \textit{outliers} - wartości odstające dalej bardzo psują nam rozk\l ady, ale zosta\l y omówione wcześniej
\item \textit{scaling} - nasze wartości są z różnych zakresów, z różną skalą, jest to niewygodne dla modelu i lepiej by to naprawi\'c. Istnieją różne sposoby, niektóre z nich są opisane na \href{https://en.wikipedia.org/wiki/Feature_scaling}{\texttt{Wikipedii}}. Są to standaryzacja, skalowanie w zależności od średniej lub min-max i podobne. Nie b\c edziemy tego robi\'c tylko dlatego, żeby nie straci\'c same dane, a raczej żeby trzyma\'c ich w początkowym stanie, bo trzymając parametry skalowania zawsze moglibyśmy ich przywróci\'c.
\item \textit{tail-heavy distributions} - jest to nasz najwi\c ekszy problem, nawet gorszy od outlier'ów, bo ich zawsze da si\c e usuną\'c, a tutaj nic nie możemy zrobi\'c, bo to jest wp\l yw na dane wprost, co powoduje zmiany ich w\l aściwości. Jednym ciekawym rozwiązaniem jest tranformacja rozk\l adu. Jest to też troch\c e opisane na \href{https://en.wikipedia.org/wiki/Data_transformation_(statistics)}{\texttt{Wikipedii}}.
\end{itemize}
\section{Analiza danych}
\quad Otóż, doszliśmy do najważniejszej i najtrudniejszej cz\c eści - analiza samych danych. Jest ona trudn ze wzgl\c edu na bardzo dużą iloś\'c regu\l i różnorodnych testów, które są potrzebne do przeprowadzenia porządnej analizy. Tym my si\c e w\l aśnie i zajmiemy.
\subsection{Wst\c ep}
\quad Musimy sobie dobrze określi\'c regu\l y, które b\c edziemy stosowa\'c dalej, żeby zawsze mie\'c możliwoś\'c odniesienia do nich.
\noindent
\quad Źród\l em regu\l\ zosta\l\ ten artyku\l\ naukowy: \textit{J. Uttley - Power Analysis, Sample Size, and Assessment of Statistical Assumptions—Improving the Evidential Value of Lighting Research} (2019). Można znaleź\'c ca\l oś\'c \href{https://www.tandfonline.com/doi/pdf/10.1080/15502724.2018.1533851}{\texttt{tutaj}}. Najbardziej interesują nas punkty 3.1 - 3.4, które są na stronach 5 - 11 w dokumencie, który znajduje si\c e pod przedstawionym linkiem. Opisują one jakie za\l ożenia stosujemy do parametrycznych testów i nieparametrycznych, i jak je sprawdzi\'c wraz z możliwymi problemami przy ich stosowaniu. Wi\c ec, ja spróbuj\c e tu krótko określi\'c jak b\c edziemy post\c epowali dalej, bazując si\c e na tym artykule naukowym.
\newpage
\noindent
\quad Odnosząc si\c e do Tabeli nr 3 w podanym artykule, możemy sobie rozdzieli\'c najcz\c eściej stosowane testy na parametryczne i nieparametryczne.
\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth]{paramtests.png}
\end{figure}
\noindent
\quad Do tej tabeli później b\c edziemy si\c e odwo\l ywa\'c, wi\c ec jest bardzo ważna dla nas.
\noindent
\quad Nie określiliśmy różnic\c e pomi\c edzy testami parametrycznymi i nieparametrycznymi. Otóż, parametryczne testy są bardziej precyzyjne, ale wymagają spe\l nienia pewnych za\l ożeń stosujących si\c e samych danych. Bez spe\l nienia ich, nie możemy stwierdzi\'c, że wynik otrzymany z tego testu jest odpowiadający rzeczywistości (nigdy tego nie możemy stwierdzi\'c, bo statystyka nigdy nie określa czegoś na 100\%, ale wtedy niepewnoś\'c jest bardzo duża). W\l aśnie możliwe problemy związane z niedotrzymaniem tych za\l ożeń są przedstawione w punkcie 3.4 w podanym artykule. Także w punkcie 3.1 jest podana informacja, że z 50 sprawdzonych artyku\l ów naukowych tylko 22\% mieli sprawdzone te za\l ożenia, co sprawia, że wyniki w nich podane mogą silnie różni\'c si\c e od rzeczywistości. To potwierdza ważnoś\'c testowania za\l ożeń i jednocześnie pokazuje jak rzadko to jest robione w akademic\c e, i to jest duży problem. W punkcie 3.5 są przedstawione przyk\l ady takiego niedotrzymania regu\l. Ale nie o tym tutaj mówimy, wi\c e b\c edziemy stara\'c si\c e w miar\c e możliwości przeprowadza\'c wszystkie możliwe badanie przed stwierdzeniem czegokolwiek.
\noindent
\quad Z tabeli nr 2 podanego artyku\l u określamy najcz\c estsze za\l ożenia testu parametrycznego.
\begin{figure}[h!]
\centering
\includegraphics[width=\textwidth]{assumptions.png}
\end{figure}
\noindent
\quad Z podanej tablicy najbardziej nas ciekawią ostatni dwa punkty, ponieważ pierwsze dwa muszą by\'c uwzgl\c ednione jeszcze przy tworzeniu samych danych, co do ostatnich dwóch, to w\l aśnie ich możemy sprawdzi\'c przed stosowaniem testów i artyku\l\ nam podaje ca\l ą list\c e sposobów na to, które b\c edziemy stosowa\'c.
\noindent
\quad Otóż, przejdźmy już do analizy podstawowej naszego zbioru danych.
<<echo=FALSE>>=
# cleaning all the variables created before
rm(list=ls())
library(corrplot)
library(purrr)
library(ggplot2)
library(hrbrthemes)
library(raster)
library(tidyverse)
library(viridis)
library(forcats)
library(dplyr)
library(ggpubr)
library(psych)
library(moments)
library(reshape)
library(reshape2)
library(MASS)
library(pastecs)
# for future visualization
density.plot.properties <- list(
geom_density(fill="darkorange", color="darkorange", alpha=1),
theme_ipsum_rc(grid='XY', base_family = 'sans'),
theme(
legend.position="none",
axis.title.x = element_blank(),
axis.title.y = element_blank()
))
plot.properties <- list(
theme_ipsum_rc(base_family = 'sans'),
theme(
legend.position="none",
axis.title.x = element_blank(),
axis.title.y = element_blank()
))
housing <- read.csv('fe_housing.csv', header=TRUE)
@
\subsection{Analiza podstawowa}
\quad Troch\c e już wcześniej zosta\l y przedstawione te dane, wi\c ec bez zb\c ednego opisu obliczmy sobie wszystkie momenty i inne elementy statystyki opisowej.
<<descriptive>>=
head(housing, 4)
summary(housing)
# all the decsriptive values (look at the skew and kurt, should see tail-heavy)
describeBy(housing)
# calculating all the moments
all.moments(housing)
# calculating all the central moments
all.moments(housing, central=TRUE)
@
\noindent
\quad \textit{Ma\l y komentarz: [1,] odpowiada za zerowy moment, [2,] za pierwszy i td...}
\noindent
\quad Jak i zosta\l o zaznaczono w komentarzu, jeżeli popatrzy\'c na wspó\l czynnik skośności i kurtoz\c e, to można zobaczy\'c, że są one w wi\c ekszości swojej dodatnie, co powoduje, że mają duży ogon, czyli są rozciągni\c ete bardziej w prawo, bardzo podobnie do rozk\l adu ex-Gaussian, czyli eksponencjalnie-normalnego. Niżej jest on przedstawiony.
\begin{figure}[h!]
\centering
\includegraphics[width=0.94\textwidth]{exgaussian.png}
\end{figure}
\noindent
\quad Narysujmy sobie jeszcze boxploty do każdej cechy.
\begin{figure}[h!]
\centering
<<descrboxplots, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
# melting the data for convenient boxplotting
melted <- melt.data.frame(housing)
melted %>% ggplot(aes(x=variable, y=value,
fill=variable)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha=0.6) +
geom_jitter(color="darkorange", size=0.05, alpha=0.1) +
theme_ipsum(base_family = 'sans') +
facet_wrap(~variable,scales="free",ncol=5, nrow=2) +
theme(
legend.position="none",
plot.title = element_text(size=20),
axis.text.x=element_text(size=12),
axis.text.y=element_text(size=12),
axis.title.x=element_blank(),
axis.title.y=element_blank()
) +
labs(title="Boxplots for each feature")
@
\includegraphics[width=0.94\textwidth]{report-descrboxplots}
\end{figure}
\noindent
\quad Widzimy dużą iloś\'c outlier'ów, które psują nam skal\c e, a także wartości różnych statystyk, wi\c e stwórzmy sobie zbiór bez nich.
<<replacingoutliers>>=
# replacing outliers with NAs to remove them while melting
outlier.replace <- function(x){
quantiles <- quantile( x, c(.003, .997 ) )
x[ x < quantiles[1] ] <- NA
x[ x > quantiles[2] ] <- NA
return(x)
}
no_outs <- housing
for(i in names(no_outs)){
no_outs[[i]] <- outlier.replace(no_outs[[i]])
}
# outliers per column
colSums(is.na(no_outs))
@
\noindent
\quad Także warto spojrze\'c na \textit{trimmed}, ponieważ to pokazuje jak zmianiają si\c e wartości po usuni\c eciu outlier'ów, i widzimy, że w przypadku niektórych cech to jest doś\'c silna zmiana.
\noindent
\quad Narysujmy teraz sobie nowe boxploty.
\begin{figure}[h!]
\centering
<<descrboxplotsnoouts, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
# melting the 'no outliers' data
melted_outs <- melt.data.frame(no_outs, na.rm=TRUE)
melted_outs %>% ggplot(aes(x=variable, y=value,
fill=variable)) +
geom_boxplot() +
scale_fill_viridis(discrete = TRUE, alpha=0.6) +
geom_jitter(color="darkorange", size=0.05, alpha=0.1) +
theme_ipsum(base_family = 'sans') +
facet_wrap(~variable,scales="free",ncol=5, nrow=2) +
theme(
legend.position="none",
plot.title = element_text(size=20),
axis.text.x=element_text(size=12),
axis.text.y=element_text(size=12),
axis.title.x=element_blank(),
axis.title.y=element_blank()
) +
labs(title="Boxplots for each feature")
@
\includegraphics[width=0.8\textwidth]{report-descrboxplotsnoouts}
\end{figure}
\noindent
\quad Narysujmy sobie jeszcze g\c estości rozk\l adów \textit{density plots}, czyli tak naprawd\c e rozmazane histogramy.
\begin{figure}[h!]
\centering
<<descrdensitplots, echo=FALSE, fig=TRUE, include=FALSE, width = 12, height = 9>>=
# density plots below
fs <- theme(axis.text.x = element_text(size=8),
axis.text.y = element_text(size=8))
p1 <- no_outs %>% ggplot(aes(x=longitude)) + density.plot.properties + fs
p2 <- no_outs %>% ggplot(aes(x=latitude)) + density.plot.properties + fs
p3 <- no_outs %>% ggplot(aes(x=housing_median_age)) + density.plot.properties + fs
p4 <- no_outs %>% ggplot(aes(x=households)) + density.plot.properties + fs
p5 <- no_outs %>% ggplot(aes(x=median_income)) + density.plot.properties + fs
p6 <- no_outs %>% ggplot(aes(x=median_house_value)) + density.plot.properties + fs
p7 <- no_outs %>% ggplot(aes(x=ocean_proximity)) + density.plot.properties + fs
p8 <- no_outs %>% ggplot(aes(x=bedrooms_per_rooms)) + density.plot.properties + fs
p9 <- no_outs %>% ggplot(aes(x=rooms_per_person)) + density.plot.properties + fs
p10 <- no_outs %>% ggplot(aes(x=people_per_household)) + density.plot.properties + fs
# drawing it all at one page
ggarrange(p1, p2, p3, p4, p5, p6, p7, p8, p9, p10,
ncol=2, nrow=5, font.label=list(size = 12, family='sans'),
labels=names(housing))
@
\includegraphics[width=\textwidth]{report-descrdensitplots}
\end{figure}
\subsection{Analiza rodzaju rozk\l adu badanych cech}
\quad Wydzielimy sobie 3 najciekawsze dla nas cechy: \textbf{bedrooms\_per\_rooms}, \textbf{median\_income} i \textbf{median\_house\_value}. Oceńmy na ile one odpowiadają rozk\l adu normalnemu, żeby potem móc przeprowadza\'c różnorodne testy na nich i mie\'c możliwoś\'c stosowania estymatorów.
\noindent
\quad Zróbmy sobie plan testowania normalności:
\begin{itemize}
\item Wizualne testowanie
\begin{itemize}
\item Histogram
\item QQ-plot
\item Boxplot
\end{itemize}
\item Oszacowanie statystyk opisowych
\begin{itemize}
\item Wspó\l czynnik skośności
\item Kurtoza
\end{itemize}
\item Statystyczne testy odchylenia od rozk\l adu nornalnego (niżej przedstawione możliwe testy)
\begin {itemize}
\item Shapiro-Wilk test
\item Kolmogorov-Smirnov test
\item Anderson-Darling test
\item D'Agostino-Pearson omnibus test
\item Jarque-Bera test
\item i inne...
\end{itemize}
\end{itemize}
\noindent
\quad Także warto zauważy\'c, że w ogólności stosujemy \textit{density plots}, ponieważ są ciekawsze i przyjemniejsze wizualnie. Ale w tym rozdziale b\c edziemy stosowa\'c histogram.
\noindent
\quad Parametr histogramu, tak zwany \textit{binwidth} może by\'c wyznaczony z regu\l y Freedman-Diaconis, która określa go jako:
\[ binwidth = \frac{2IQR}{\sqrt[3]{n}}\]
\noindent
\quad Także b\c edziemy stosowa\'c QQ-plot. Niżej jest wygląd tego wykresu dla danych o rozk\l adzie normalnym (czyli to jest prosta).
\begin{figure}[h!]
\centering
<<qqplotnorm, echo=TRUE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
ggplot(data.frame(d = rnorm(1000)), aes(sample=d)) + stat_qq()
@
\includegraphics[width=\textwidth]{report-qqplotnorm}
\end{figure}
\noindent
\quad Podczas stosowania statystyki opisowej do sprawdzenia normalności b\c edziemy pos\l ugiwa\'c si\c e $z$-statystyką bazującą si\c e na wspó\l czynniku skośności i kurtozie. Jest ona wyznaczona jako:
\[ z = \frac{skew}{2SE} \quad z = \frac{kurtosis}{2SE} \]
\noindent
\quad I w końcu, ze wszystkich pokazanych wyżej możliwych testów na normalnoś\'c b\c edziemy stosowali \textit{Shapir-Wilk test}, ponieważ jest on bardziej czutliwy i mocny, zgodnie z pracami cytowanymi w artykule. A także warto zauważy\'c, że powienien by\'c przeprowadzany na próbkach o rozmiarze nie wi\c ekszym od 50. Jest to powiązane, że w opisie tego testu od 1965 roku przez Shapiro i Wilk'a byli stosowane takie ilości, bo nie mieli dost\c epu do wi\c ekszych. Jest to też opisane w artykule. O ile mamy dużą iloś\'c rekordów, to nie możemy wprost polega\'c na przypadkowości wybierania próbki. Dlatego b\c edziemy robi\'c po 3 testy dla różnych próbek 50-elementowych.
\noindent
\quad Podczas przeprowadzenia testu zak\l adamy, że \textbf{$H_0$} - hipoteza zerowa - stwierdza, że rozk\l ad jest normalnym, a \textbf{$H_1$} - hipoteza alternatywna - stwierdza, że rozk\l ad takim nie jest.
\subsubsection{\textit{bedrooms\_per\_rooms}}
\begin{figure}[h!]
\centering
<<histbpr, echo=TRUE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
gap <- IQR(no_outs$bedrooms_per_rooms, na.rm=TRUE) /
(length(na.omit(no_outs$bedrooms_per_rooms)) ^ (1/3))
no_outs %>% ggplot(aes(x=bedrooms_per_rooms)) + plot.properties +
geom_histogram(binwidth = gap) +
labs(title='Bedrooms per rooms distribution')
@
\includegraphics[width=0.99\textwidth]{report-histbpr}
\end{figure}
\noindent
\quad Z histogramu wida\'c ten ogon, i możemy już przewidywa\'c, że rozk\l ad nie b\c edzie rozk\l adem normalnym, ale przeprowadźmy dalszą analiz\c e.
\begin{figure}[h!]
\centering
<<qqplotbpr, echo=FALSE, fig=TRUE, include=FALSE, width = 8, height = 6>>=
no_outs %>% ggplot(aes(sample=bedrooms_per_rooms))+stat_qq() +
labs(title='Bedrooms per rooms QQ-plot')
@
\includegraphics[width=\textwidth]{report-qqplotbpr}
\end{figure}
\noindent
\quad Widzimy bardzo duże odchylenie od prostej, wi\c ec możemy by\'c pewnymi, że to nie jest rozk\l ad normalny, ale przeprowadźmy dalszą analiz\c e, żeby sobie dodatkowo zobaczy\'c jak b\c edzie si\c e zachowywa\l\ rozk\l ad nie normalny. Dla innych cech potem, jeżeli zobaczymy, że zbiór nie odpowiada rozk\l adowi normalnemu, to nie b\c edziemy kontynuowali obliczeń dla niego.
\noindent
\quad Nast\c epnym punktem jest boxplot, ale zosta\l\ już wyżej przedstawiony. Punkty pomi\c edzy 1 i 3 kwartylem są doś\'c podobnie roz\l ożone odnośnie mediany, wąsy są doś\'c blisko d\l ugości $1.5IQR$, chociaż wyraźnie wida\'c, że górny jest troszeczk\c e wi\c ekszy. Co pokazuje nienormalnoś\'c - to iloś\'c outlier'ów i ich tendencja pojawia\'c si\c e z góry, czyli nie ma dolnych outrlier'ów.
\noindent
\quad Dalej testujemy wartoś\'c $z$-statystyki bazującej si\c e na wspó\l czynniku skośności i kurtozie. Powinien by\'c testowany zgodnie z artyku\l em na próbkach nie wi\c ekszych od 200.
<<normkurtbpr>>=
stat.desc(sample(no_outs$bedrooms_per_rooms, 150), norm=TRUE)
@