Domača naloga pri predmetu MTRD
2023-11-04
Ime in priimek: Denis Marinšek
1. Izbor podatkov, vir podatkov, cilj analize ipd.
library(psych)
podatki <- force(sat.act) #Uvozim podatke
head(podatki)## gender education age ACT SATV SATQ
## 29442 2 3 19 24 500 500
## 29457 2 3 23 35 600 500
## 29498 2 3 20 21 480 470
## 29503 1 4 27 26 550 520
## 29504 1 2 33 31 600 550
## 29518 1 5 26 28 640 640colnames(podatki) <- c("Spol", "Izobrazba", "Starost", "Uspeh", "Verbal", "Kvant") #Preimenujem spremenljivke
head(podatki) #Prikažem prvih 6 vrstic tabele s podatki## Spol Izobrazba Starost Uspeh Verbal Kvant
## 29442 2 3 19 24 500 500
## 29457 2 3 23 35 600 500
## 29498 2 3 20 21 480 470
## 29503 1 4 27 26 550 520
## 29504 1 2 33 31 600 550
## 29518 1 5 26 28 640 640Opis spremenljivk:
- Spol: 1:M, 2:Z
- Izobrazba: Dosežena izobrazba (kategorije)
- Uspeh: Dosežen splošen uspeh
- Verbal: Rezultat na preizkusu iz verbalnih sposobnosti
- Kvant: Rezultat na preizkusu iz kvantitativnih znanj
Vir podatkov: RStudio, knjižnica Psych, podatki: sat.act.
Cilj analize: Podatke sem izbral, ker bi rad analiziral razlike v splošnem uspehu, verbalnih ter kvantitativnih znanjih med spoloma.
2. Urejanje podatkov, izvedba opisne statistike in njena razlaga
podatki$SpolF <- factor(podatki$Spol, #Ustvarim faktor
levels = c(1, 2),
labels = c("M", "Z"))
library(tidyr)
podatki <- drop_na(podatki) #Odstranim menjakjoče vrednosti
podatki$ID <- seq(1, nrow(podatki)) #Ustvarim spremenljivko ID
summary(podatki[ , c(-1, -2, -8)]) #Prikažem opisno statistiko## Starost Uspeh Verbal Kvant SpolF
## Min. :13.00 Min. : 3.00 Min. :200.0 Min. :200.0 M:245
## 1st Qu.:19.00 1st Qu.:25.00 1st Qu.:550.0 1st Qu.:530.0 Z:442
## Median :22.00 Median :29.00 Median :620.0 Median :620.0
## Mean :25.64 Mean :28.55 Mean :612.3 Mean :610.2
## 3rd Qu.:29.00 3rd Qu.:32.00 3rd Qu.:700.0 3rd Qu.:700.0
## Max. :65.00 Max. :36.00 Max. :800.0 Max. :800.0Razlaga: Ocena povprečne starosti za osebe, zbrane v vzorec, znaša 25,64 let. V vzorcu imamo 442 žensk.
#Prikažem porazdelitev verbalnih sposobnosti s histogramom
hist(podatki$Verbal,
main = "Porazdelitev rezultatov iz verbalnih sposobnosti",
xlab = "Dosežen rezutat iz verbalnih sposobnosti",
ylab = "Frekvenca")
Razlaga: Ugotavljam, da je porazdelitev asimetrična v levo.
3. Raziskovalno vprašanje, na katerega bom odgovoril s pomočjo statističnega preverjanja domnev: Ali se dosežena rezultata pri Kvantiativnih znanjih med spoloma razlikujeta?
#S pomočjo histogramov ločeno prikažem porazdelitev kvantitativnih znanj po spolu
library(ggplot2)
Moski <- ggplot(podatki[podatki$SpolF == "M", ], aes(x = Kvant)) +
theme_linedraw() +
geom_histogram(binwidth = 10) +
ylab("Frekvenca") +
ggtitle("Moski")
Zenske <- ggplot(podatki[podatki$SpolF == "Z", ], aes(x = Kvant)) +
theme_linedraw() +
geom_histogram(binwidth = 10) +
ylab("Frekvenca") +
ggtitle("Zenske")
library(ggpubr)
ggarrange(Moski, Zenske,
ncol = 2, nrow = 1)
# Izvedem Shapiro-Wilkov preizkus za preverbo normalnosti
library(dplyr)
library(rstatix)
podatki %>%
group_by(SpolF) %>%
shapiro_test(Kvant)## # A tibble: 2 × 4
## SpolF variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 M Kvant 0.943 0.0000000335
## 2 Z Kvant 0.971 0.000000101Razlaga: Ker se porazdelitvi razlikujeta od pričakovane normalne porazdelitve (tako grafično kot na podlagi Shapiro-Wilkovih preizkusov, obrakrat p < 0,001), izberem neparametrični preizkus.
library(psych)
describeBy(podatki$Kvant, podatki$SpolF) #Opisna statistika po skupinah, da bom primerjal mediane##
## Descriptive statistics by group
## group: M
## vars n mean sd median trimmed mad min max range skew kurtosis
## X1 1 245 635.87 116.02 660 645.53 94.89 300 800 500 -0.72 -0.12
## se
## X1 7.41
## ------------------------------------------------------------
## group: Z
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 442 596 113.07 600 602.21 133.43 200 800 600 -0.58 0.13 5.38#Wilcoxnov preizkus vsote rangov
wilcox.test(podatki$Kvant ~ podatki$SpolF,
paired = FALSE,
correct = FALSE,
exact = FALSE,
alternative = "two.sided")##
## Wilcoxon rank sum test
##
## data: podatki$Kvant by podatki$SpolF
## W = 65921, p-value = 2.223e-06
## alternative hypothesis: true location shift is not equal to 0H0: Lokaciji porazdelitev uspeha pri kvantiativnih znanjih sta enaki
H1: Lokaciji porazdelitev uspeha pri kvantiativnih znanjih se razlikujeta
Razlaga: Na podlagi vzorčnih podatkov zavrnemo ničelno domnevo pri p < 0,001. Ugotavljamo, da se lokaciji porazdelitev uspeha pri kvantiativnih znanjih razlikujeta. Glede na vrednosti median ugotavljamo, da so moški dosegli boljše rezultate kot ženske.
library(effectsize)#Izračunam velikost učinka
effectsize(wilcox.test(podatki$Kvant ~ podatki$SpolF,
paired = FALSE,
correct = FALSE,
exact = FALSE,
alternative = "two.sided"))## r (rank biserial) | 95% CI
## --------------------------------
## 0.22 | [0.13, 0.30]interpret_rank_biserial(0.22)## [1] "medium"
## (Rules: funder2019)Razlaga: Našli smo srednje velike razlike v doseženem rezultatu pri kvantitativnih znanjih.
4. Raziskovalno vprašanje, na katerega bom odgovoril s pomočjo regresijske analize: Ali lahko na podlagi rezultatov testov predvidimo spol osebe?
#Logistična regresija
fit <- glm(SpolF ~ Uspeh + Verbal + Kvant,
family = binomial,
data = podatki)library(car)
vif(fit) #Preverba multikolinearnosti## Uspeh Verbal Kvant
## 1.753972 1.919772 2.049930Razlaga: V modelu je prisotna določena stopnja multikolinearnost, vendar bomo vse spremenljivke upoštevali v končnem modelu, saj so vse vrednosti še vedno nižje od 5.
podatki$StdOstanki <- rstandard(fit) #Izračun std. ostankov
podatki$Cook <- cooks.distance(fit) #Izračun Cookovih razdalj
hist(podatki$StdOstanki,
main = "Histogram standardiziranih ostankov",
ylab = "Frekvenca",
xlab = "Standardizirani ostanki")
Razlaga: Pričakovana oblika porazdelitve, opazimo tudi da ni osamelcev.
head(podatki[order(-podatki$Cook), c("ID", "Cook")]) ## ID Cook
## 279 279 0.02450469
## 434 434 0.01782366
## 182 182 0.01710302
## 544 544 0.01571290
## 498 498 0.01476430
## 99 99 0.01293209library(dplyr)
podatki <- podatki %>%
filter(!ID == 279) #Odstranim enoto ID 279Razlaga: Odstranimo enoto z oznako 279, saj predstavlja enoto z visokim vplivom.
fit0 <- glm(SpolF ~ 1,
family = binomial,
data = podatki)
fit <- glm(SpolF ~ Uspeh + Verbal + Kvant,
family = binomial,
data = podatki)anova(fit0, fit, test = "Chi") #Primerjava dveh logističnih modelov## Analysis of Deviance Table
##
## Model 1: SpolF ~ 1
## Model 2: SpolF ~ Uspeh + Verbal + Kvant
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 685 893.03
## 2 682 861.81 3 31.223 7.63e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Razlaga: Model z upoštevanimi tremi pojasnjevalnimi spremenljivkami se bolje prilega podatkom kot osnovni model brez pojasnjevalnih spremenljivk (p < 0,001).
summary(fit) #Povzetek modela##
## Call:
## glm(formula = SpolF ~ Uspeh + Verbal + Kvant, family = binomial,
## data = podatki)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.757182 0.552877 3.178 0.00148 **
## Uspeh 0.028645 0.022629 1.266 0.20557
## Verbal 0.002435 0.001021 2.386 0.01705 *
## Kvant -0.005634 0.001077 -5.232 1.68e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 893.03 on 685 degrees of freedom
## Residual deviance: 861.81 on 682 degrees of freedom
## AIC: 869.81
##
## Number of Fisher Scoring iterations: 4exp(cbind(RO = fit$coefficients, confint.default(fit))) #Razmerja obetov## RO 2.5 % 97.5 %
## (Intercept) 5.7960823 1.9612115 17.1294992
## Uspeh 1.0290589 0.9844158 1.0757266
## Verbal 1.0024382 1.0004346 1.0044457
## Kvant 0.9943823 0.9922859 0.9964831Razlaga: Če se rezultat pri verbalnih spretnostih poveča za 1 točko, je obet, da je oseba ženskega spola, 1,0026 kratnik začetnega obeta (p = 0,012), ob predpostavki, da so ostale pojasnjevalne spremenljvike nespremenjene. [razložite celoten model]
podatki$OceneVerjet <- fitted(fit) #Izračun ocen verjetnosti# Izračun psevdo R2 za ugotovitev odstotka pravilno uvrščenih enot.
podatki$Uvrstitev <- ifelse(test = podatki$OceneVerjet < 0.50,
yes = "NE",
no = "DA")
podatki$UvrstitevF <- factor(podatki$Uvrstitev,
levels = c("NE", "DA"),
labels = c("NE", "DA"))
razvrst_tabela <- table(podatki$SpolF, podatki$UvrstitevF)
razvrst_tabela##
## NE DA
## M 36 208
## Z 20 422Psevdo_R2_fit <- (razvrst_tabela[1, 1] + razvrst_tabela[2, 2] )/ nrow(podatki)
Psevdo_R2_fit## [1] 0.6676385Razlaga: Z modelom pravilno uvrstimo 66 % oseb.