Prva domača naloga

Djenita Mustafoska

Prikaz pripravljenih podatkov za analizo:

library(NHANES)
nh_data <- NHANES
podatki <- nh_data[c("Gender", "Age", "Weight", "Height", "Poverty")]
head(podatki)

##   Gender Age Weight Height Poverty
## 1   male  34   87.4  164.7    1.36
## 2   male  34   87.4  164.7    1.36
## 3   male  34   87.4  164.7    1.36
## 4   male   4   17.0  105.4    1.07
## 5 female  49   86.7  168.4    1.91
## 6   male   9   29.8  133.1    1.84

Prikazani podatki so iz knjižnice NHANES (National Health and Nutrition Examination Survey). Za našo nalogo so bili izbrani samo nekateri podatki (spremenljivke) katere se je zdelo smiselno analizirati.

Spremenljivke:

Gender je kategorična spremenljivka, ki predstavlja spol posameznika udeleženega v raziskavi, ta je moški ali ženski, označeno v tabeli z (male, female).
Age je številska spremenljivka, ki predstavlja starost posameznika udeleženega v raziskavi.
Weight je številska spremenljivka, ki predstavlja težo posameznika, zaokrožena na eno decimalko.
Height je številska spremenljivka, ki predstavlja višino posameznika, zaokrožena na eno decimalko.
Poverty je kazalnik revščine gospodinjstva, izračunan kot skupni družinski dohodek deljen s pragom revščine, ki ustreza velikosti in sestavi tega gospodinjstva (indeks manjši od 1 pomeni, da je dohodek gospodinjstva pod pragom revščine in velja za revno) in je tudi številska spremenljivka.
Za te podatke je potrebno naložiti paket “NHANES”, to naredimo z ukazom install.packages(“NHANES”).

Cilj analize:

Naredili bomo preprosto analizo, kjer bomo poskušali ugotoviti kateri spol je v povprečju višji in tudi težji (saj sklepamo, da sta ti dve lastnosti korelirani). In s tem da imamo še podatke o kazalniku revščine gospodinjstva za posameznika, bomo preverili kateri spol je v povprečju premožnejši.

Čiščenje podatkov:

Imamo zelo velik vzorec 10000 opazovanj (prevelik, da bi naredili nekatere teste, kot so Shapirov Wilkov test za preizkus normalnosti), torej naši analizi ne bo škodilo če nekatere podatke odstranimo, še posebej tiste, ki nam ne povejo veliko. Odstranimo vsa opazovanja, ki imajo v kateri koli spremenljivki vrednost “NA”. Spremenimo tudi spremenljivko spola v faktor, kjer bo po novem 0 predstavljala moški spol in 1 ženski spol.

podatki <- na.omit(podatki)
library(dplyr)
podatki <- podatki %>%
  mutate(Gender = factor(Gender, levels = c("male", "female"), labels = c(0, 1)))
summary(podatki)

##  Gender        Age           Weight           Height         Poverty     
##  0:4461   Min.   : 2.0   Min.   :  9.60   Min.   : 83.6   Min.   :0.000  
##  1:4475   1st Qu.:19.0   1st Qu.: 57.90   1st Qu.:156.9   1st Qu.:1.260  
##           Median :37.0   Median : 73.80   Median :166.1   Median :2.720  
##           Mean   :37.6   Mean   : 72.94   Mean   :161.9   Mean   :2.818  
##           3rd Qu.:54.0   3rd Qu.: 89.70   3rd Qu.:174.7   3rd Qu.:4.760  
##           Max.   :80.0   Max.   :230.70   Max.   :200.4   Max.   :5.000

Po čiščenju podatkov imajo naše spremenljivke 8936 opazovanj.

Vizualizacija podatkov

par(mfrow = c(2, 2))

x = range(80,205)

x_2 = range(0,205)

y = range(0, 600)

hist(podatki$Height[podatki$Gender == 0], 
     main = "Histogram moške višine",
     xlab = "Višina",
     ylab = "Frekvenca",
     col = "lightblue",
     border = "black",
     breaks = 50,
     xlim = x,
     ylim = y)

hist(podatki$Weight[podatki$Gender == 0], 
     main = "Histogram moške teže",
     xlab = "teža",
     ylab = "Frekvenca",
     col = "lightblue",
     border = "black",
     breaks = 50,
     xlim = x_2,
     ylim = y)

hist(podatki$Height[podatki$Gender == 1], 
     main = "Histogram ženske višine",
     xlab = "Višina",
     ylab = "Frekvenca",
     col = "pink",
     border = "black",
     breaks = 50,
     xlim = x,
     ylim = y)

hist(podatki$Weight[podatki$Gender == 1], 
     main = "Histogram ženske teže",
     xlab = "teža",
     ylab = "Frekvenca",
     col = "pink",
     border = "black",
     breaks = 50,
     xlim = x_2,
     ylim = y)

Analiza podatkov

Za začetek preverimo koreliranost med višino in težo, da vidimo smiselnost našega vzorca, saj te dve spremenljivki očitno morata biti korelirani (odvisni) med sabo in nekako pričakujemo, da zavrnemo ničelno hipotezo, da spremenljivki nista korelirani (sta neodvisni). Gremo zapisati to hipotezo in narediti korelacijski test za to hipotezo.

H₀ : Med dvema spremenljivkama ni statistično pomembne linearne korelacije (koeficient korelacije je enak 0).
H₁ : Med dvema spremenljivkama obstaja statistično pomembna linearna korelacija (koeficient korelacije ni enak 0).

cor.test(podatki$Height,podatki$Weight, method = "pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  podatki$Height and podatki$Weight
## t = 107.55, df = 8934, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7419506 0.7600271
## sample estimates:
##       cor 
## 0.7511296

Dobimo p-vrednost = 2.2e-16 < 0.05, kar pomeni, da ničelno hipotezo zavrnemo, in ne moremo trditi, da spremenljivki nista korelirani. Tudi p-vrednost je zelo majhna, kar nakazuje zelo močno povezanost med tema dvema spremenljivkama, s čimer smo zadovoljni.

library(ggplot2)

data <- data.frame(podatki$Weight, podatki$Height)

# Create a stacked bar chart
ggplot(data, aes(x = podatki$Weight, y = podatki$Height, color = podatki$Gender)) +
  geom_point(size = 0.3, shape = 19) +
  labs(title = "Razsevni diagram teže in višine, moških in žensk pri 8936 opazovanjih.",
       x = "teža",
       y = "višina") +
  theme_minimal() +
  scale_color_manual(values = c("blue", "red"), labels = c("moški", "ženske")) +  
  theme(legend.title = element_blank())

Tudi razsevni diagram nakazuje, da sta spremenljivki odvisni.

Histogram nakazuje, da so moški višji in težji od žensk, zanima nas če je to res.

Prvo s Shapiro-Wilkovim testom preverimo če so višine in teže porazdeljene normalno:

H₀ : Vzorec je porazdeljen normalno.
H₁ : Vzorec ni porazdeljen normalno.

shapiro.test(podatki$Height[podatki$Gender == 0])

## 
##  Shapiro-Wilk normality test
## 
## data:  podatki$Height[podatki$Gender == 0]
## W = 0.75299, p-value < 2.2e-16

shapiro.test(podatki$Height[podatki$Gender == 1])

## 
##  Shapiro-Wilk normality test
## 
## data:  podatki$Height[podatki$Gender == 1]
## W = 0.76829, p-value < 2.2e-16

shapiro.test(podatki$Weight[podatki$Gender == 0])

## 
##  Shapiro-Wilk normality test
## 
## data:  podatki$Weight[podatki$Gender == 0]
## W = 0.96525, p-value < 2.2e-16

shapiro.test(podatki$Weight[podatki$Gender == 1])

## 
##  Shapiro-Wilk normality test
## 
## data:  podatki$Weight[podatki$Gender == 1]
## W = 0.97399, p-value < 2.2e-16

Pri vseh preizkusih Shapiro-Wilkovega testa zavrnemo H₀, da ima vzorec normalno porazdelitev. To pomeni, da t-test o preizkušanju vrednosti aritmetične sredine med dvema vzorcema ni primeren (saj za ta test mora biti vzorec porazdeljen normalno).

Naredimo Wilcoxonov test za neodvisna vzorca. Preizkušamo naslednjo domnevo:

H₀ : Lokacija moške višine/teže vzorca je večja ali enaka lokaciji ženske višine/teže vzorca.
H₁ : Lokacija moške višine/teže vzorca je manjša lokaciji ženske višine/teže vzorca

wilcox.test(podatki$Height ~ podatki$Gender,
       paired = FALSE,
       correct = FALSE, exact = FALSE,
       alternative = "less")

## 
##  Wilcoxon rank sum test
## 
## data:  podatki$Height by podatki$Gender
## W = 15633110, p-value = 1
## alternative hypothesis: true location shift is less than 0

wilcox.test(podatki$Weight ~ podatki$Gender,
       paired = FALSE,
       correct = FALSE, exact = FALSE,
       alternative = "less")

## 
##  Wilcoxon rank sum test
## 
## data:  podatki$Weight by podatki$Gender
## W = 12689973, p-value = 1
## alternative hypothesis: true location shift is less than 0

V obeh primerih je p-vrednosti = 1, torej ne zavrnemo H₀, kar se nam tudi zdi smiselno glede na naš vzorec.

Možno je, da nam anomalije (mlade osebe ali osebe z motnjo v razvoju) pokvarijo vzorec saj ima naš vzorec dolg lev rep in zgleda približno normalno porazdeljen. Opravimo še t-test za naslednjo hipotezo:

H₀ : Aritmetična sredina vzorca višine/teže moških je večja od aritmetične sredine vzorca višine/teže žensk.
H₁ :Aritmetična sredina vzorca višine/teže moških je manjša od aritmetične sredine vzorca višine/teže žensk.

t.test(Height ~ Gender, data = podatki, alternative = "less")

## 
##  Welch Two Sample t-test
## 
## data:  Height by Gender
## t = 25.174, df = 8346.7, p-value = 1
## alternative hypothesis: true difference in means between group 0 and group 1 is less than 0
## 95 percent confidence interval:
##      -Inf 11.14503
## sample estimates:
## mean in group 0 mean in group 1 
##        167.1666        156.7052

t.test(Weight ~ Gender, data = podatki, alternative = "less")

## 
##  Welch Two Sample t-test
## 
## data:  Weight by Gender
## t = 17.381, df = 8746, p-value = 1
## alternative hypothesis: true difference in means between group 0 and group 1 is less than 0
## 95 percent confidence interval:
##     -Inf 11.0041
## sample estimates:
## mean in group 0 mean in group 1 
##        77.97671        67.92407

Dobimo p-vrednost = 1 >0.05, torej tako kot prej tudi pri tem testu ne zavrnemo H₀ hipoteze. Lahko sklepamo, da so moški v povprečju višji in težji od žensk.

Probajmo ugotoviti kateri spol je bolj premožen v povprečju.

Postavimo domnevo za Wilcoxonov test:

H₀ : Ni statistično pomembne razlike v porazdedlitvi indikatorju revščine med spoloma.
H₁ : So statistično pomembne razlike v porazdedlitvi indikatorju revščine med spoloma.

wilcox.test(Poverty ~ Gender, data = podatki, paired = FALSE,
       correct = FALSE, exact = FALSE, alternative = "two.sided")

## 
##  Wilcoxon rank sum test
## 
## data:  Poverty by Gender
## W = 10293575, p-value = 0.01003
## alternative hypothesis: true location shift is not equal to 0

Oboje stranski preizkus nam da p-vrednost = 0.015515 < 0.05, kar pomeni da zavrnemo H₀ hipotezo in sprejmemo alternativno torej so statistično pomembne razlike med spoloma.

Preverimo še če ženske v povprečju služijo manj kot moški. Postavimo hipotezo:

H₀ : Lokaciji porazdelitve indikatorja revščine sta enaki.
H₁ : Lokaciji porazdelitve indikatorja revščine nista enaki.

wilcox.test(Poverty ~ Gender, data = podatki, alternative = "less")

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  Poverty by Gender
## W = 10293575, p-value = 0.995
## alternative hypothesis: true location shift is less than 0

Dobimo p-vrednost = 0.995 > 0.05 in ne zavrnemo ničelno hipotezo. Predpostavimo, da moški v popvprečju živijo v bolj premožnem gospodinjstvu kot ženske. Ta podatek (indikator revščine) si lahko nekateri osebki delijo, zato mogoče ni najprimernejše kazalo kateri spol je premožnejši, saj si lahko osebki različnega spola delijo gospodinjstvo mi pa nevemo točno kako je bila ta anketa izvedena.

Poglejmo si še porazdelitev bogatstva med spoli.

par(mfrow = c(1, 2))

hist(podatki$Poverty[podatki$Gender == 0], 
     main = "Histogram ind. revščine za moške",
     xlab = "indikator revščine",
     ylab = "Frekvenca",
     col = "chartreuse3",
     border = "black",
     breaks = 50,
     xlim = range(0,5),
     ylim = range(0,1500))

hist(podatki$Poverty[podatki$Gender == 1], 
     main = "Histogram ind. revščine za ženske",
     xlab = "indikator revščine",
     ylab = "Frekvenca",
     col = "khaki2",
     border = "black",
     breaks = 50,
     xlim = range(0,5),
     ylim = range(0,1500))

library(ggplot2)

ggplot(podatki, aes(x = Gender, y = Poverty, fill = as.factor(Gender))) +
  geom_boxplot() +
  labs(x = "spol", y = "indeks revščine") +
  ggtitle("graf kvantilov za indikator revščine med spoloma") +
  scale_fill_manual(values = c("lightblue", "pink"), labels = c("moški", "ženske")) +
  theme(legend.title = element_blank())

Vizualni podatki malo nakazujejo neenakost v indikatorju revščine med spoloma, kar smo tudi v naši raziskavi ugotovili.