Domača naloga s programom R pri predmetu MTRD
Urban Celin
2024-01-08
1. RAZISKOVALNO VPRAŠANJE: ALI SE INDEKS TELESNE MASE V ZDA GLEDE NA REGIJO/OBMOČJE RAZLIKUJE?
1.1 PODATKI
SPREMENLJIVKE:
- Starost: dopolnjena starost v letih
- Spol: biološki spol (male = M, female = Z)
- BMI: vrednost indeksa telesne mase v kilogramih na kvadratni meter (Kg/m2). Normalen razpon se giblje od 18.5 do 24. 9
- StOtrok: Število otrok, ki jih zavzema zdravstveno zavarovanje
- Kadilec: Status kadilca (yes = DA, no = NE)
- Regija: Lokacijsko območje v ZDA (southeast = JV, southwest = JZ, northeast = SV, northwest = SZ)
- Stroski: Znesek računa enkratne zdravstvene obravnave izražen v ameriških dolarjih ($)
VIRI PODATKOV : Choi, M. (2017). Medical Cost Personal Datasets. Pridobljeno 23. 12. 2023 s https://www.kaggle.com/datasets/mirichoi0218/insurance
UVOZ IN UREJANJE SEKUNDARNIH PODATKOV:
podatki <- read.table("./insurance.csv",
header = TRUE,
sep = ",",
dec = ".") #Uvozim podatke v Coma-separated values obliki
colnames(podatki) <- c("Starost", "Spol", "BMI", "StOtrok", "Kadilec", "Regija", "Stroski")
#Spremenljivkam določim imena
podatki$ID <- seq(1, nrow(podatki)) #Dodam ID-je posameznim enotam
head(podatki) #Izpišem prvih šest vrednosti ## Starost Spol BMI StOtrok Kadilec Regija Stroski ID
## 1 19 female 27.900 0 yes southwest 16884.924 1
## 2 18 male 33.770 1 no southeast 1725.552 2
## 3 28 male 33.000 3 no southeast 4449.462 3
## 4 33 male 22.705 0 no northwest 21984.471 4
## 5 32 male 28.880 0 no northwest 3866.855 5
## 6 31 female 25.740 0 no southeast 3756.622 6podatki$SpolF <- factor(podatki$Spol, #Spremenljivko Spol pretvorim v faktor
levels = c( "male", "female"),
labels = c("M", "Z"));
podatki$KadilecF <- factor(podatki$Kadilec, #Spremenljivko Kadilec pretvorim v faktor
levels = c( "yes", "no"),
labels = c("DA", "NE"));
podatki$RegijaF <- factor(podatki$Regija, #Spremenljivko regija pretvorim v faktor
levels = c("southeast", "southwest", "northeast", "northwest"),
labels = c("JV", "JZ", "SV", "SZ"));
str(podatki) #Preverim strukturo posameznih spremenljivk ## 'data.frame': 1338 obs. of 11 variables:
## $ Starost : int 19 18 28 33 32 31 46 37 37 60 ...
## $ Spol : chr "female" "male" "male" "male" ...
## $ BMI : num 27.9 33.8 33 22.7 28.9 ...
## $ StOtrok : int 0 1 3 0 0 0 1 3 2 0 ...
## $ Kadilec : chr "yes" "no" "no" "no" ...
## $ Regija : chr "southwest" "southeast" "southeast" "northwest" ...
## $ Stroski : num 16885 1726 4449 21984 3867 ...
## $ ID : int 1 2 3 4 5 6 7 8 9 10 ...
## $ SpolF : Factor w/ 2 levels "M","Z": 2 1 1 1 1 2 2 2 1 2 ...
## $ KadilecF: Factor w/ 2 levels "DA","NE": 1 2 2 2 2 2 2 2 2 2 ...
## $ RegijaF : Factor w/ 4 levels "JV","JZ","SV",..: 2 1 1 4 4 1 1 4 3 4 ...set.seed(1700) #določim seme, da vedno delam na istih podatkih
podatki <- podatki[sample(nrow(podatki), 150), ] #Ker gre za obširen dataset umetno zmanjšam vzorec na n = 150OPIS PODATKOV:
summary(podatki[colnames(podatki) %in% c("Starost", "SpolF", "BMI", "StOtrok", "KadilecF", "RegijaF", "Stroski")]) ## Starost BMI StOtrok Stroski SpolF KadilecF
## Min. :18.0 Min. :17.20 Min. :0.00 Min. : 1136 M:75 DA: 35
## 1st Qu.:27.0 1st Qu.:25.80 1st Qu.:0.00 1st Qu.: 5016 Z:75 NE:115
## Median :40.5 Median :30.30 Median :1.00 Median :10097
## Mean :40.0 Mean :30.33 Mean :1.12 Mean :14190
## 3rd Qu.:53.0 3rd Qu.:34.27 3rd Qu.:2.00 3rd Qu.:19359
## Max. :64.0 Max. :45.54 Max. :4.00 Max. :48173
## RegijaF
## JV:40
## JZ:40
## SV:32
## SZ:38
##
## #izpišem opisno stat. s pomočjo funkcije summary, da lahko med drugim preverim tudi frekvence v faktorjih RAZLAGA ZA SPREM. STAROST IN FAKTOR SPOL: Najmlajša oseba zajeta v vzorec je bila stara 18 let, najstarejša pa 64 let. Prvi kvantil predstavlja 27 let, tretji kvantil pa 53 let. Mediana starosti je znašala 40,5 let - natanko polovica enot je mlajša ali stara natanko 40,5 let, ostale enote so starejše od 40,5 let. Aritmetična sredina znaša natanko 40 let. V vzorec je zajetih 75 moških in 75 žensk.
library(psych)
describe(podatki [, c(-2, -5, -6, -8, -9, -10, -11)]) #izpisem opisno statistko s pomočjo funkcije describe,## vars n mean sd median trimmed mad min max
## Starost 1 150 40.00 14.12 40.50 39.97 20.02 18.0 64.00
## BMI 2 150 30.33 5.93 30.30 30.39 6.41 17.2 45.54
## StOtrok 3 150 1.12 1.09 1.00 1.01 1.48 0.0 4.00
## Stroski 4 150 14190.36 12430.02 10096.53 12071.05 8624.74 1136.4 48173.36
## range skew kurtosis se
## Starost 46.00 -0.01 -1.26 1.15
## BMI 28.34 -0.06 -0.56 0.48
## StOtrok 4.00 0.57 -0.75 0.09
## Stroski 47036.96 1.32 0.75 1014.91#smiselno iz izpisa odstranim faktorje in ID-je RAZLAGA ZA SPREMENLJIVKO STROŠKI: V vzorec je zajetih 150 oseb, aritmetična sredina enkratnih stroškov zdravstvene obravnave znaša 14190,36 $, standardni odklon znaša 12430, 02 $, mediana znaša 10096, 53 $, aritmetična sredina brez ekstremnih vrednosti znaša 12071, 05 $, mediana absolutne deviacije znaša 8624, 74 $, najnižja vrednost stroškov znaša 1136, 40 $, najvišja vrednost stroškov znaša 48173, 36 $, variacijski razmik znaša 47036, 96 $, porazdelitev stroškov zdravstvene obravnave je asimetrična v desno ter pretežno koničasta, standardna napaka znaša 1014, 91 $.
BodyMI <- hist(podatki$Stroski, #Prikažem histogram porazdelitve stroškov zdravstvenih storitev
main = "Porazdelitev stroškov zdravstvene oskrbe ",
xlab = "Stroški zdravstvene oskrbe ($) ",
ylab = "Frekvenca",
freq=TRUE,
breaks = 20,
col="darkmagenta", #polnilo
border="black") #obroba 
RAZLAGA: Glede na grafični prikaz lahko sklepam,
da je porazdelitev stroškov asimetrična v desno ter pretežno koničasta.
Gre za pričakovano obliko porazdelitve.
1.2 ANALIZA PODATKOV
CILJI ANALIZE: Na podlagi sekundarnih vzorčnih podatkov analizirati ali se indeks telesne teže v ZDA razlikuje med posameznimi regijami/območji.
library(psych)
describeBy(podatki$BMI, podatki$RegijaF) #describe opisna statistika, grupirano po posameznih regijah ##
## Descriptive statistics by group
## group: JV
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 40 32.36 5.93 33.33 32.41 6.6 21.12 45.54 24.42 -0.06 -0.8 0.94
## ------------------------------------------------------------
## group: JZ
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 40 31.74 5.78 33.4 31.93 6.08 17.4 42.9 25.5 -0.34 -0.48 0.91
## ------------------------------------------------------------
## group: SV
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 32 28.03 6.21 27.55 27.84 6.83 17.2 39.71 22.52 0.19 -1.03 1.1
## ------------------------------------------------------------
## group: SZ
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 38 28.65 4.81 28.88 28.76 4.23 17.39 37.72 20.33 -0.32 -0.2 0.78RAZLAGA: Na podlagi primerjave centralnih tendenc med posameznimi področji ugotavljam, da lahko pričakujem statistično značilne razlike med posameznimi področji. Iz opisanega sledi, da je razlike smiselno statistično preveriti. Vsebinsko pričakujem, da se bo od drugih regij razlikovala predvsem jugovzhodna regija, saj so zvezdni državi Florida in Georgia, ki so del le te eni izmed vodilnih držav v povezavi z prekomerno telesno težo.
#install.packages('dplyr')
#install.packages('rstatix')
library(rstatix)
library(dplyr)
podatki %>%
group_by(RegijaF) %>% #Izvedem test normalnosti BMI-ja na posameznem območju s pomočjo Shapiro-Wilkovega testa normalnosti
shapiro_test(BMI)## # A tibble: 4 × 4
## RegijaF variable statistic p
## <fct> <chr> <dbl> <dbl>
## 1 JV BMI 0.976 0.551
## 2 JZ BMI 0.977 0.588
## 3 SV BMI 0.967 0.409
## 4 SZ BMI 0.976 0.574RAZLAGA: Shapirov test normalnost: H0 = podatki sledijo normalni porazdelitvi, H1 = podatki ne sledijo normalni porazdelitvi. Na podlagi Shapiro-Wilkovega testa sklepam, da se spremenljivka BMI v vseh štirih področjih porazdeljuje normalno (p > 0,05 v vseh štirih skupinah).
#install.packages("tidyverse")
#install.packages('ggplot2')
library(ggplot2) #aktiviram knjižnico ggplot2
JugoV <- ggplot(podatki[podatki$RegijaF == "JV", ], aes(x = BMI)) + #Jugovzhodno področje
theme_linedraw() +
geom_histogram(binwidth = 3, color="#F4DF4EFF",
fill="#949398FF", linewidth = 1.2) +
geom_density(alpha = 0.5, fill="#FF6655") +
xlab("Indeks telesne mase")+
ylab("Frekvenca") +
ggtitle("Porazdelitev BMI-ja na jugovzhodu")
JugoZ <- ggplot(podatki[podatki$RegijaF == "JZ", ], aes(x = BMI)) + #Jugozadhodno področje
theme_linedraw() +
geom_histogram(binwidth = 3, color="#F4DF4EFF",
fill="#949398FF", linewidth = 1.2) +
geom_density(alpha = 0.5, fill="#FF6655") +
xlab("Indeks telesne mase")+
ylab("Frekvenca") +
ggtitle("Porazdelitev BMI-ja na jugozahodu")
SeveroV <- ggplot(podatki[podatki$RegijaF == "SV", ], aes(x = BMI)) + #Severovzhodna regija
theme_linedraw() +
geom_histogram(binwidth = 3, color="#E69A8DFF",
fill="#5F4B8BFF", linewidth = 1.2) +
geom_density(alpha = 0.5, fill="#FF6655") +
xlab("Indeks telesne mase")+
ylab("Frekvenca") +
ggtitle("Porazdelitev BMI-ja na severovzhodu")
SeveroZ <- ggplot(podatki[podatki$RegijaF == "SZ", ], aes(x = BMI)) + #Severozahodna regija
theme_linedraw() +
geom_histogram(binwidth = 3, color="#E69A8DFF",
fill="#5F4B8BFF", linewidth = 1.2) +
geom_density(alpha = 0.5, fill="#FF6655") +
xlab("Indeks telesne mase")+
ylab("Frekvenca") +
ggtitle("Porazdelitev BMI-ja na severozahodu")
#install.packages('ggpubr')
library(ggpubr)
ggarrange(JugoV, JugoZ, SeveroV, SeveroZ, #grafične prikaze postavim v dve vrsti in dva stolpca
ncol = 2, nrow = 2)
RAZLAGA: Na podlagi grafičnega prikaza sklepam, da porazdlitve BMI-ja vizualno izgledajo dokaj normalne na vseh štirih področjih.Vizualno je še najbolj problematična porazdelitev severovzhodnega predela. Kljub temu Shapiro-Wilkov test potrjuje normalnost v vseh skupinah zato smatram, da je predpostavka o normalni porazdlitvi obravnavane spremenljivke znotraj raziskovanih skupin izpolnjena.
#install.packages('car')
library(car)
leveneTest(podatki$BMI, podatki$RegijaF) #Izvedem Levene-ov test homogenosti varianc glede na štiri področja ## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.3915 0.2478
## 146RAZLAGA: Levene-ov test homogenosti varianc: H0 = variance so homogene, H1 = variance niso homogene. Na podlagi rezultata testa zaključim, da je predpostavka o homoskedastičnosti izpolnjena (F = 1.3915, p > 0.05)
rezultatiAOV <- aov(BMI ~ RegijaF, #na vzorčnih podatkov izvedem parametrični test analize varianc ANOVA
data = podatki)
summary (rezultatiAOV) #izpis povzetka ## Df Sum Sq Mean Sq F value Pr(>F)
## RegijaF 3 521 173.76 5.367 0.00156 **
## Residuals 146 4727 32.38
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#install.packages('effectsize')
library(effectsize) #Preverim velikost učinka
eta_squared(rezultatiAOV)## # Effect Size for ANOVA
##
## Parameter | Eta2 | 95% CI
## -------------------------------
## RegijaF | 0.10 | [0.03, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].interpret_eta_squared(0.10, rules = "cohen1992") #interpretacija velikost učinka## [1] "small"
## (Rules: cohen1992)RAZLAGA: Na podlagi parametričnega testa analize varianc ugotavljam, da med področji/regijami v ZDA pri indeksu telesne mase prihaja do statistično značilnih razlik (vsaj ena izmed primerjanih aritmetičnih sredin se razlikuje od ostalih) (p < 0,002). Razlike so vsebinsko majhne (Eta2 = 0.10).
library(rstatix)
pairwise_t_test(BMI ~ RegijaF, #Izvedem Post hoc test, da ugotovim med katerimi skupinami prihaja do razlik
paired = FALSE, #gre za skupine
p.adjust.method ="bonferroni", #uporabim Bonferronijev popravek, da se izognem napihnjenim p-jem
data = podatki)## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 p p.signif p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <chr> <dbl> <chr>
## 1 BMI JV JZ 40 40 0.621 ns 1 ns
## 2 BMI JV SV 40 32 0.00163 ** 0.00975 **
## 3 BMI JZ SV 40 32 0.00682 ** 0.0409 *
## 4 BMI JV SZ 40 38 0.00452 ** 0.0271 *
## 5 BMI JZ SZ 40 38 0.0179 * 0.107 ns
## 6 BMI SV SZ 32 38 0.652 ns 1 nsRAZLAGA: Na podlagi Post-Hoc testov ugotavljam, da do statistično značilnih razlik med indeksom telesne mase prihaja med jugovzhodnim predelom ter severnovzhodnim predelom (p < 0,01 ), med jugozahodnim ter severozahodnim predelom (p < 0,005) ter med jugovzhodnim in severozahodnim predelom (p < 0,03). Razlike med ostalimi tremi pari so statistično neznačilne.
kruskal.test(BMI ~ RegijaF, #Opravim alternativni neparametrični Kruskal-Wallisov test
data = podatki)##
## Kruskal-Wallis rank sum test
##
## data: BMI by RegijaF
## Kruskal-Wallis chi-squared = 14.145, df = 3, p-value = 0.002715library(effectsize)
kruskal_effsize(BMI ~ RegijaF, #Ponovno ocenim velikost učinka
data = podatki)## # A tibble: 1 × 5
## .y. n effsize method magnitude
## * <chr> <int> <dbl> <chr> <ord>
## 1 BMI 150 0.0763 eta2[H] moderateRAZLAGA: Na podlagi Kruskal-Wallisovega testa lahko trdim, da med BMI-jem med posameznimi regijami prihaja do statistično značilnih razlik (vsaj ena preučevana lokacija porazdelitve se razlikuje od preostalih) (Chi squared = 14,145, p < 0.003 ). Razlike so vsebinsko gledano srednje velike (Eta2 = 0.0763).
library(rstatix)
wilcox_test(BMI ~ RegijaF, #Post hoc test, da lahko opredelim razlike med posameznimi skupinami
paired = FALSE, #gre za skupine
p.adjust.method = "bonferroni", #Bonferronijev popravek, da se izognem napaki prve vrste
data = podatki)## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <chr>
## 1 BMI JV JZ 40 40 840. 0.707 1 ns
## 2 BMI JV SV 40 32 878 0.007 0.043 *
## 3 BMI JV SZ 40 38 1038 0.006 0.033 *
## 4 BMI JZ SV 40 32 868 0.01 0.06 ns
## 5 BMI JZ SZ 40 38 1007 0.014 0.082 ns
## 6 BMI SV SZ 32 38 560. 0.579 1 nsRAZLAGA: Na podlagi Wilcox Post-Hoc testa ugotavljam, da do statistično značilnih razlik pri indeksih telesne mase prihaja med jugovzhodnim in severovzdhodnim predelom ZDA (p < 0.05) ter med jugovzhodnim predelom in severozahodnim predelom ZDA (p < 0.04).
1.3 UGOTOVITVE
Ker predpostavki o normalnosti obravnavane spremenljivke znotraj
skupin ter o homogenosti varianc niso bile kršene, obravnavana
spremenljivka pa je bila številska in ni vsebovala osamelcev, kot bolj
relevantne rezultate vidim tiste, ki jih ponujajo parametrični testi
(saj imajo le ti v temu primeru večjo statistično moč). Na podlagi
vzorčnih podatkov sem s pomočjo analize varianc (ANOVA) potrdil, da med
skupinami (predeli ZDA) prihaja do statistično značilnih razlik (vsaj
ena aritmetična sredina se je razlikovala od ostalih) (p < 0,02, F =
5.367). Velikost učinka je bila majhna (Eta2 = 0.10). Nadaljno sem na
podlagi vzorčnih podatkov uspel s pomočjo Post-Hoc parametričnih testov
statistično potrdil, da do statistično značilnih razlik pri indeksu
telesne mase prihaja med jugovzhodnim predelom ter severnovzhodnim
predelom ZDA (p < 0,01 ), med jugozahodnim ter severozahodnim
predelom ZDA (p < 0,005) ter med jugovzhodnim in severozahodnim
predelom ZDA (p < 0,03). Razlike med ostalimi tremi pari so
statistično neznačilne. S slednjim sem uspešno odgovoril na zastavljeno
raziskovalno vprašanje - razlike v indeksu telesne mase so med
nekaterimi področji v ZDA prisotne. Prav tako so ugotovitve v skladu z
pričakovanji.
2. RAZISKOVALNO VPRAŠANJE: KATERI DEJAVNIKI VPLIVAJO NA VIŠINO STROŠKOV ENKRATNE ZDRAVSTVENE OBRAVNAVE V ZDA?
2.1 PODATKI
SPREMENLJIVKE:
- Starost: dopolnjena starost v letih
- Spol: biološki spol (male = M, female = Z)
- BMI: vrednost indeksa telesne mase v kilogramih na kvadratni meter (Kg/m2). Normalen razpon se giblje od 18,5 do 24, 9
- StOtrok: Število otrok, ki jih zavzemna zdravstveno zavarovanje
- Kadilec: Status kadilca (yes = DA, no = NE)
- Regija: Lokacijsko območje v ZDA (southeast = JV, southwest = JZ, northeast = SV, northwest = SZ)
- Stroski: Znesek računa enkratno opravljene zdravstvene obravnave izražen v ameriških dolarjih ($)
VIRI PODATKOV : Choi, M. (2017). Medical Cost Personal Datasets. Pridobljeno 23. 12. 2023 s https://www.kaggle.com/datasets/mirichoi0218/insurance
UVOZ IN UREJANJE SEKUNDARNIH PODATKOV: Ponoven uvoz ter urejanje sekundarnih podatkov ni potreben. Statistična analiza bo temeljila na že predhodno uvoženih, urejenih in obdelanih podatkih. Po potrebi podatki sproti smiselno urejam tokom analize.
CILJI ANALIZE: Preveriti ali lahko s pomočjo
izbranih dejavnikov predvidimo (pojasnimo) stroške enkratne zdravstvene
oskrbe v ZDA. Pričakujem, da bosta porast indeksa telesne mase in
starosti negativno vplivala na stroške (jih v povprečju višala). Prav
tako pričakujem, da bodo ne-kadilci v primerjavi z kadilci imeli v
povprečju manjše stroške. Med bio. spoli ne pričakujem razlik. Na
podlagi prej izvedene analize bo predvidoma prihajalo do razlik tudi med
posameznimi območji v ZDA.
2.2 ANALIZA PODATKOV
fit0 <- lm(Stroski ~ 1, #linearna regresija
data = podatki)
fit1 <- lm(Stroski ~ Starost, #Blokovna postopna gradnja modela s pomočjo dodajanje pojas. sprem.
data = podatki)
fit2 <- lm(Stroski ~ Starost + BMI,
data = podatki)
fit4 <- lm(Stroski ~ Starost + SpolF + BMI + StOtrok + KadilecF + RegijaF,
data = podatki)
anova(fit0, fit1, fit2, fit4, test = "Chi") #Primerjava modelov ## Analysis of Variance Table
##
## Model 1: Stroski ~ 1
## Model 2: Stroski ~ Starost
## Model 3: Stroski ~ Starost + BMI
## Model 4: Stroski ~ Starost + SpolF + BMI + StOtrok + KadilecF + RegijaF
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 149 2.3021e+10
## 2 148 2.0332e+10 1 2.6889e+09 < 2.2e-16 ***
## 3 147 1.9091e+10 1 1.2417e+09 1.172e-08 ***
## 4 141 5.3818e+09 6 1.3709e+10 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1RAZLAGA: Na podlagi analize variance (primerjave modelov) ugotavljam, da se model fit4 naboljše prilega podatkom (p < 0,0001).
podatki$StOstanki <- rstandard(fit4) #standardizirani ostanki
podatki$CooksD <- cooks.distance(fit4) #Cookove razdalje
hist(podatki$StOstanki, # histrogram stand. ostankov
xlab = "Standardizirani ostanki",
ylab = "Frekvenca",
main = "Histogram standardiziranih ostankov",
col="orange",
border="black")
hist(podatki$CooksD, #histogram Cookovih razdalj
xlab = "Cookove razdalje",
ylab = "Frekvenca",
main = "Histogram Cookovih razdalj", #Glavni naslov
col="beige", # polnilo
border="black") #obroba
RAZLAGA: Na podlagi grafičnega prikaza ugotavljam, da sta v podatkih prisotna vsaj dva osamelca ter vsaj tri enoti z visokim vplivom. Porazdelitvi sta pričakovani.
head(podatki[order(-podatki$StOstanki), c(8,12) ], 8) #izpišem osem najvišjih vrednosti stand. ostankov ## ID StOstanki
## 431 431 3.396764
## 469 469 3.063688
## 103 103 2.737505
## 600 600 2.660400
## 492 492 2.157583
## 1259 1259 2.140975
## 246 246 2.091307
## 697 697 2.065155head(podatki[order(-podatki$CooksD), c(8,13) ], 8) #Izpišem osem najvišjih vrednosti Cookovih razdalj## ID CooksD
## 431 431 0.07766977
## 103 103 0.05737784
## 600 600 0.05181484
## 469 469 0.04821770
## 492 492 0.04338254
## 1259 1259 0.03929638
## 430 430 0.02794327
## 697 697 0.02614847library(dplyr)
podatki2 <- podatki %>% # filtriram podatke, da se znebim osamelcev/enot s visokim vplivom
filter(!ID %in% c(431, 103, 600, 469 ))RAZLAGA: Na podlagi izpisa ugotovim, da enota z ID-jem 431 pradstavlja tako osamelca kot enoto z visokim vplivom. Ob enem enoti z ID-jem 103 in 600 predstavljata enoti z visokim vplivom, enota z ID-jem 469 pa je samo osamelec. Omenjene enote iz vzorca izključim (n = 146)
library(car)
scatterplot(Stroski ~ Starost + SpolF , #Prikažem linearni vpliv Starosti na stroške glede na spol
ylab = "Stroški zdravstvene oskrbe ",
xlab = "Starost",
main = "Preverjanje lin. med Stroški in Starostjo glede na Spol ",
smooth = FALSE,
regLine = FALSE, #Brez regresijske premice
boxplot = FALSE, #Brez grafikona kvantilov
legend = FALSE, #Brez legende
data = podatki2)
podatkistev <- podatki2[, c("Stroski", "Starost", "BMI", "StOtrok")] #samo številske spre.
library(car)
scatterplotMatrix(podatkistev,
smooth = FALSE,
main = "Preverjanje linearnosti in porazdelitev spremenljivk v modelu") #naslov
quantile(podatki2$Stroski, 0.95) # Poiščem vrednosti 95. centila Stroškov ## 95%
## 41858.47podatki2$StroskiW <- ifelse (test = podatki2$Stroski > 41807.7,
yes = 41858.47, #Zamenjam vrednosti
no = podatki2$Stroski)fit5 <- lm(StroskiW ~ Starost + SpolF + BMI + KadilecF + RegijaF, #Testiram modela
data = podatki2)
fit6 <- lm(StroskiW ~ Starost + SpolF + BMI + StOtrok + KadilecF + RegijaF,
data = podatki2)
anova(fit5, fit6, test = "Chi")## Analysis of Variance Table
##
## Model 1: StroskiW ~ Starost + SpolF + BMI + KadilecF + RegijaF
## Model 2: StroskiW ~ Starost + SpolF + BMI + StOtrok + KadilecF + RegijaF
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 138 3638442511
## 2 137 3602660491 1 35782020 0.2434##install.packages('lmtest')
library(lmtest) #Iz knjižnice lmtest izvozim Rainbow test ter le tega izvedem, da presodim ali v modelu velja linearnost
raintest(fit5)##
## Rainbow test
##
## data: fit5
## Rain = 1.1201, df1 = 73, df2 = 65, p-value = 0.3215RAZLAGA: S pomočjo grafičnih prikazov preverjam normalnosti porazdelitev ter do neke mere tudi linearno povezanost spremenljivk. Ker je spremenljivka Stroski asimetrična v desno na njej izvedem winsorizacijo. Opažam, da tudi pri spremenljivki StOtrok prihaja do težav. Ker na podlagi analize varianc zaključim, da slednja v pomembni meri ne prispeva k kvaliteti prileganja modela podatkov le to iz modela odstranim. Na koncu predpostavko o linearnosti preverim tudi s pomočjo dodatnega formalnega Rainbow testa. Slednji preverja, če je dobro linearno prileganje možno doseči v podvzorcu, ki je lociran na sredini podatkov. Če je dobro prileganje dejansko možno doseči, spremenljivke med sabo niso linearne povezane. Ničelna hipoteza navaja, da med spremenljivkami v osnovnemu modelu vlada linearna povezanost (odstranitev podatkov signifikantno ne pridonese k boljšemu prileganju), alternativna pa, da je slednja povezanost znotraj modela med spremenljivkami nelinearna (odstranitev podatkov povzroči bistveno boljše prileganje). Na podlagi testa ničelne hipoteze ne morem zavrniti. Predpostavka o linearnosti ni kršena (p > 0.05).
#install.packages('olsrr')
library(olsrr)
ols_test_breusch_pagan(fit5) #Izvedem Breusch-Paganov test, da preverim predpostavko o homoskedastičnosti ##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## ------------------------------------
## Response : StroskiW
## Variables: fitted values of StroskiW
##
## Test Summary
## ------------------------------
## DF = 1
## Chi2 = 16.57334
## Prob > Chi2 = 4.68043e-05RAZLAGA: Breusch-Pagan test: H1 = variance so homogene, H2 = variance so heterogene. Na podlagi Breusch-Paganovega testa ugotavljam, da je predpostavka o homoskedastičnosti kršena (p < 0,001).
##install.packages('lmtest')
library(lmtest) #Uvozim knjižnico v kateri se med drugih nahajajo alternativni testi za preveranje homogenosti
gqtest(fit5, order.by = ~Starost + BMI, data = podatki2, fraction = 29 ) #Starost bolj varira -> prvo mesto##
## Goldfeld-Quandt test
##
## data: fit5
## GQ = 3.4126, df1 = 51, df2 = 50, p-value = 1.301e-05
## alternative hypothesis: variance increases from segment 1 to 2#Goldfeld-Quandtov test, 29 = število centralih enot, ki jih test ignorira (običajno 20% vzorca)
bptest(fit5, ~ Starost + SpolF + BMI + KadilecF + RegijaF, data = podatki2) ##
## studentized Breusch-Pagan test
##
## data: fit5
## BP = 19.522, df = 7, p-value = 0.0067#Whitov test za preverjanje homoskedastičnosti podatki2$StOstanki <- rstandard(fit5)
podatki2$StdOcenVred <- scale(fit5$fitted.values)
library(car)
scatterplot (StOstanki ~ StdOcenVred,
main = "Preverjanje predpostavke homoskedističnosti",
ylab = "Standardizirani ostanki",
xlab = "Standardizirane ocenjene vrednosti",
boxplot = FALSE,
smooth = FALSE,
data = podatki2)
RAZLAGA: Predpostavko o homoskedastičnosti preverim tudi z alternativnim Goldfeld-Quandtovim parametričnim testom homogenosti varianc. Na podlagi testa ničelno hipotezo zavračam (p < 0,001). Rezultat testa je v skladu z rezultati Breusch-Paganovega testa. Uporabim tudi Whitov test konstantnosti varianc -> h0 = variance so konstantne, h1 = variance so ne-konstantne. Tudi v temu primeru ničelno hipotezo zavrnem (p < 0,005), rezultat testa je skladen s preostalimi dognanji. Ker so vsi testi konsistentni lahko predpostavljam, da do napak pri specifikaciji spremenljivk predvidoma ni prišlo. V nasprotnem primeru je Whitov test pogosto nekonsistenten z preostalimi alternativnimi testi za preverjanje homoskedastičnosti. Nazadnje predpostavko preverim tudi grafično. Kot kaže je predpostavka o homoskedastičnosti kršena. Da zagotovim nepristranskost bom moral v regresijski model vključiti Whitove robustne standardne napake. Predhodno poskusim tudi s procesom tranformacije kljub temu, da ima le ta v osnovi drugačen namen.
shapiro.test (podatki2$StOstanki) ##
## Shapiro-Wilk normality test
##
## data: podatki2$StOstanki
## W = 0.95333, p-value = 7.909e-05RAZLAGA: Predpostavko o normalnosti standardiziranih napak preverjam s pomočjo Shapiro-Wilkovega testa normalnosti, ki ga opravim na standardiziranih ostankih: H0 = stand. ostanki se porazdeljujejo normalno, H1 = stand. ostanki se ne porazdeljujejo normalno. Ničelno hipotezo zavračam (p < 0,001) iz česar sledi, da prihaja do kršitve predpostavke. Slednje ne bi smelo biti problematično glede na velikost vzorca (n = 146).
#install.packages('car')
library(car)
vif(fit5)## GVIF Df GVIF^(1/(2*Df))
## Starost 1.042401 1 1.020981
## SpolF 1.038169 1 1.018906
## BMI 1.152555 1 1.073571
## KadilecF 1.017408 1 1.008666
## RegijaF 1.166413 3 1.025988mean(vif(fit5))## [1] 1.171004RAZLAGA: Na podlagi VIF statistike zaključim, da do težav z multikolinearnostjo med spremenljivkami, ki sem jih vključil v model ne prihaja. Posamezne VIF vrednosti so krepto pod 5, povprečna VIF statistika je blizu 1. Predpostavka je potemtakem izpolnjena
#Logaritmiram odvisno in neodvisne sprem.
fit5loglog <- lm(log1p(StroskiW) ~ log1p(Starost) + SpolF + log1p(BMI) + KadilecF + RegijaF,
data = podatki2)
fit5log <- lm(log1p(StroskiW) ~ Starost + SpolF + BMI + KadilecF + RegijaF,
data = podatki2) #Logaritmiram samo odvisno sprem.
anova (fit5log, fit5loglog, test ="Chi") #Primerjam prileganje modelov ## Analysis of Variance Table
##
## Model 1: log1p(StroskiW) ~ Starost + SpolF + BMI + KadilecF + RegijaF
## Model 2: log1p(StroskiW) ~ log1p(Starost) + SpolF + log1p(BMI) + KadilecF +
## RegijaF
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 138 20.429
## 2 138 20.011 0 0.41784library(olsrr)
ols_test_breusch_pagan(fit5log) #Še enkrat preverim homoskedastičnost ##
## Breusch Pagan Test for Heteroskedasticity
## -----------------------------------------
## Ho: the variance is constant
## Ha: the variance is not constant
##
## Data
## -------------------------------------------
## Response : log1p(StroskiW)
## Variables: fitted values of log1p(StroskiW)
##
## Test Summary
## ----------------------------
## DF = 1
## Chi2 = 0.2204863
## Prob > Chi2 = 0.6386697RAZLAGA: Ker predpostavka o homoskedastičnosti ni bila izpolnjena sem model fit5 najprej logaritmiral. V osnovi sem logaritmiral tako odvisno kot neodvisne spremenljivke. Nato sem logaritmiral samo odvisno spremenljivko ter primerjal prileganje modelov. Odločim se za model fitlog. Ponovno izvedem Breusch-Paganov test, ničelne hipoteze ne morem zavrniti (p > 0.05). Tranformacija je odpravila težave z heterogenostjo varianc.
library(estimatr)
fit_robust <- lm_robust(StroskiW ~ Starost + SpolF + BMI + KadilecF + RegijaF, #reg. s pomočjo rob. stand. napak
se_type = "HC1", #uporabim HC1 popravek
data = podatki2)
summary(fit5log)##
## Call:
## lm(formula = log1p(StroskiW) ~ Starost + SpolF + BMI + KadilecF +
## RegijaF, data = podatki2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.95946 -0.19356 -0.01271 0.12293 1.39658
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.098804 0.221039 36.640 < 2e-16 ***
## Starost 0.034986 0.002328 15.029 < 2e-16 ***
## SpolFZ 0.110276 0.064895 1.699 0.091516 .
## BMI 0.021322 0.005754 3.705 0.000305 ***
## KadilecFNE -1.450276 0.075233 -19.277 < 2e-16 ***
## RegijaFJZ -0.011756 0.087617 -0.134 0.893460
## RegijaFSV 0.143221 0.096722 1.481 0.140953
## RegijaFSZ 0.093745 0.091067 1.029 0.305089
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3848 on 138 degrees of freedom
## Multiple R-squared: 0.8198, Adjusted R-squared: 0.8107
## F-statistic: 89.69 on 7 and 138 DF, p-value: < 2.2e-16summary(fit_robust)##
## Call:
## lm_robust(formula = StroskiW ~ Starost + SpolF + BMI + KadilecF +
## RegijaF, data = podatki2, se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 2609.9 3406.87 0.7661 4.449e-01 -4126.5 9346.3 138
## Starost 299.2 33.63 8.8968 2.883e-15 232.7 365.7 138
## SpolFZ 1179.6 856.15 1.3778 1.705e-01 -513.3 2872.5 138
## BMI 492.6 89.37 5.5114 1.692e-07 315.8 669.3 138
## KadilecFNE -21821.0 1253.19 -17.4124 1.213e-36 -24299.0 -19343.1 138
## RegijaFJZ -1001.8 1132.32 -0.8847 3.778e-01 -3240.7 1237.1 138
## RegijaFSV 1337.2 1398.87 0.9559 3.408e-01 -1428.8 4103.2 138
## RegijaFSZ 174.4 1305.75 0.1336 8.939e-01 -2407.4 2756.3 138
##
## Multiple R-squared: 0.8197 , Adjusted R-squared: 0.8106
## F-statistic: 81.05 on 7 and 138 DF, p-value: < 2.2e-16RAZLAGA:S pomočjo vključevanja White-ovih
robustnih napak oblikujem tudi robustni model. Odločim se za HC1
popravek, saj smo le tega uporabljali tokom predmeta. Sam menim, da je v
primeru heteroskedastičnosti manj pristranski HC3 popravek. Na koncu
izpišem povzetek tako logaritmiranega modela (fit5log) kot robustnega
modela (fit_robust). Odločim se za interpretacijo robustnega
modela.
2.3 UGOTOVITVE
S pomočjo multiple linearne regresije sem prišel do relevantnih ugotovitev. Ugotavljam, da na višino stroškov zdravstvene obravnave statistično značilno vplivajo tri dejavniki. Interpretacije parcialnih koeficientov so naslednje:
Če se oseba postara za 1 leto, se v povprečju strošek zdravstvene obravnave poveča za 299, 2 $ ob predpostavki, da ostale neodvisne spremenljivke ostanejo nespremenjene (p < 0,0001).
Če se indeks telesne mase (BMI) osebe poveča za eno indeksno točko, se v povprečju strošek zdravstvenih obravnave poveča za 492, 6 $ ob predpostavki, da ostale neodvisne spremenljivke ostanejo nespremenjene (p < 0,0001).
Osebe, ki niso kadilci imajo ob danih vrednostih ostalih neodvisnih spremenljivk v povprečju za 21821 $ manj stroškov pri zdravstvenih obravnavah kot osebe, ki so kadilci (p < 0,0001).
Ostalih parcilanih koeficientov ne smem interpretirati, ker le ti niso statistično značilni. Populacijski koeficient determinacije R squared znaša 0.8197. Iz slednjega lahko sklepamo, da smo z uporabljenim modelom uspeli pojasniti 81, 97 % variacije odvisne spremenljivke (stroškov zdravstvene obravnave). Na podlagi F statistike (h0 = vsi parcialni koeficienti so enaki nič, h1 = vsaj en parcialni koeficient ni enak nič) zavrnem ničelno hipotezo (p < 0.0001, F = 81.05). Regresijski model je statistično značilen. Rezultati so v skladu z pričakovanji - med spoloma ni signifikantnih razlik, tako višji BMI kot višja starost negativni vplivajo na stroške zdravniške obravnave, med ne-kadilci in kadilci so prisotne razlike. Kljub drugačnim pričakovanjem razlike med regijami niso signifikantne.