1 Cel projektu

Wykorzystanie danych pacjentów do oszacowania średnich wydatków na opiekę medyczną przy użyciu regresji dla przewidzianych segmentów populacji.

1.1 Zbiór danych

W niniejszym projekcie korzystamy ze zbioru danych insurance.csv. Zestaw ten został zainspirowany książką Machine Learning with R autorstwa Brett Lantz. Dane zawierają 1338 przykładów beneficjentów aktualnie zarejestrowanych w planie ubezpieczeniowym, z charakterystykami wskazującymi cechy pacjenta, a także całkowitymi wydatkami medycznymi obciążającymi plan ubezpieczeniowy na rok kalendarzowy. Wśród danych medycznych i genetycznych, możemy znaleźć także informacje o obszarze mieszkalnym beneficjenta (w USA) oraz czy jest palaczem tytoniu.

Źródło: kaggle.com

Zmienne

Tabela 1. Opis zmiennych
Zmienna	Opis
age	wiek (lata)
sex	płeć
bmi	wskaźnik masy ciała (\(kg/m^{2})\)
children	liczba dzieci objętych ubezpieczeniem zdrowotnym / liczba osób pozostających na utrzymaniu
smoker	czy beneficjent jest palaczem
region	obszar mieszkalny beneficjenta w USA: Northeast; Southeast; Southwest; Northwest
charges	indywidualne koszty leczenia rozliczane przez ubezpieczenie zdrowotne (USD)

1.2 Wczytanie danych

Tabela 2. Przedstawienie danych
age	sex	bmi	children	smoker	region	charges
19	female	27.900	0	yes	southwest	16884.924
18	male	33.770	1	no	southeast	1725.552
28	male	33.000	3	no	southeast	4449.462
33	male	22.705	0	no	northwest	21984.471
32	male	28.880	0	no	northwest	3866.855
31	female	25.740	0	no	southeast	3756.622

2 Czyszczenie danych

Zaczniemy od przedstawienia struktury danych.

Tabela 3. Struktura danych
Zmienna	Typ zmiennej	Wartości
age	integer	19, 18, 28, 33, 32, 31
sex	character	female, male, male, male, male, female
bmi	double	27.9, 33.77, 33, 22.705, 28.88, 25.74
children	integer	0, 1, 3, 0, 0, 0
smoker	character	yes, no, no, no, no, no
region	character	southwest, southeast, southeast, northwest, northwest, southeast
charges	double	16884.924, 1725.5523, 4449.462, 21984.47061, 3866.8552, 3756.6216

Zbiór składa się z 1338 obserwacji i posiada 7 kolumn.

Skontrolujmy, także czy nie występują wartości NA.

##      age      sex      bmi children   smoker   region  charges 
##    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE    FALSE

Po wstępnej analizie uznaliśmy, że zmienne sex, smoker i region powinny być zmiennymi kategorycznymi. Przekształciliśmy je do odpowiedniej postaci, przy pomocy funkcji factor. Podzieliśmy także zmienne na dwa podzbiory, testowy oraz uczący się, co będzie kluczowe przy analizie modelowej.

set.seed(35)
split <- sample.split(insurance, SplitRatio = 2/3)
insurance_train <- subset(insurance, split == TRUE)
insurance_test <- subset(insurance, split == FALSE)

3 Eksploracyjna analiza danych

3.1 Podstawowe statystyki

Ponieważ program nie zawiera funkcji, która oblicza modę stworzyliśmy taką, która ją tworzy:

W pierwszej kolejności przedstawiamy częstotliwość występowania nienumerycznych zmiennych:

Tabela 4. Liczebność zmiennych nienumerycznych

sex	Freq
female	380
male	384

smoker	Freq
no	608
yes	156

children	Freq
0	326
1	182
2	129
3	100
4	16
5	11

region	Freq
northeast	186
northwest	188
southeast	209
southwest	181

Następnie przedstawiamy podstawowe statystyki dla zmiennych numerycznych, których wartości przedstawione są na poniższej tabeli:

Age

Charges

BMI

Tabela 5. Podstawowe statystyki

	Female	Male
Min	18.00	18.00
Max	64.00	64.00
Q1	27.00	27.00
Median	41.00	40.00
Q3	52.00	52.00
Mode	19.00	18.00
Mean	39.94	39.78
Sd	14.16	14.39
IQR	25.00	25.00
Sx	12.50	12.50
Var %	0.35	0.36
IQR Var %	0.61	0.62
Skewness	-0.01	0.01
Kurtosis	-1.28	-1.27

	a	b
Min	18.00	18.00
Max	64.00	64.00
Q1	27.00	27.00
Median	41.00	39.50
Q3	52.00	51.00
Mode	18.00	19.00
Mean	39.95	39.50
Sd	14.27	14.29
IQR	25.00	24.00
Sx	12.50	12.00
Var %	0.36	0.36
IQR Var %	0.61	0.61
Skewness	-0.02	0.09
Kurtosis	-1.30	-1.19

	Female	Male
Min	1615.77	1135.94
Max	63770.43	62592.87
Q1	5000.98	4759.51
Median	9741.51	9756.70
Q3	14489.39	17437.11
Mode	3756.62	1725.55
Mean	12742.27	14115.86
Sd	11038.17	13302.09
IQR	9488.41	12677.60
Sx	4744.20	6338.80
Var %	0.87	0.94
IQR Var %	0.97	1.30
Skewness	1.70	1.43
Kurtosis	2.70	1.10

	a	b
Min	1135.94	12829.46
Max	36580.28	63770.43
Q1	4238.65	20931.48
Median	7727.05	33904.10
Q3	11411.18	42111.81
Mode	1725.55	27808.73
Mean	8536.59	32514.81
Sd	5835.74	11983.24
IQR	7172.53	21180.33
Sx	3586.26	10590.17
Var %	0.68	0.37
IQR Var %	0.93	0.62
Skewness	1.52	0.14
Kurtosis	3.45	-1.04

	Female	Male
Min	16.82	15.96
Max	48.07	53.13
Q1	26.40	27.08
Median	30.20	30.80
Q3	34.58	35.44
Mode	29.92	30.88
Mean	30.59	31.29
Sd	6.13	6.24
IQR	8.18	8.36
Sx	4.09	4.18
Var %	0.20	0.20
IQR Var%	0.27	0.27
Skewness	0.29	0.33
Kurtosis	-0.26	0.18

	a	b
Min	15.96	17.20
Max	53.13	52.58
Q1	26.78	26.11
Median	30.50	30.20
Q3	34.78	35.58
Mode	32.11	28.31
Mean	30.94	30.95
Sd	6.08	6.62
IQR	8.00	9.47
Sx	4.00	4.74
Var %	0.20	0.21
IQR Var%	0.26	0.31
Skewness	0.26	0.45
Kurtosis	-0.04	-0.05

^a No smoker

^b Smoker

Na podstawie powyższych tabel można przedstawić pierwsze spostrzeżenia analizy wzrokowej, które zostaną później sprawdzone poprzez pogłębione anliazy:

age, możemy stwierdzić, że średni wiek pacjenta, to 39.94 lata (kobiety) i 39.78 lat (mężczyźni), najstarszy beneficjent, korzystający z planu ubezpieczeniowego miał 64 lata, a najmłodszy 18 lat. 50% pacjentów mieściło się w przedziale 27 - 52 lat. Odchylenie standardowe wynosi 14.16 lat (kobiety) i 14.39 lat (mężczyźni).
sex, podział ze względu na płeć jest w miarę zrównoważony, tj. 380 kobiet i 384 mężczyzn.
bmi, średni przyrost masy ciała dla mężczyzn wyniósł 31.29 \(kg/m^{2}\) (otyłość I stopnia), a dla kobiet 30.59 \(kg/m^{2}\)., maksymalna wartość tego wskaźnika to 53.13 \(kg/m^{2}\) (otyłość III stopnia), a minimalna 15.96 (wygłodzenie). 50% pacjentów płci męskiej mieściło się w przedziale 27.08 - 35.44 \(kg/m^{2}\) (nadwaga/otyłość I stopnia), a płci żeńskiej w przedziale 26.40 - 34.58 \(kg/m^{2}\). Odchylenie standardowe wynosi odpowiednio 6.13 i 6.24 kg/m^2 dla kobiet i mężczyzn.
children, najczęściej żadne z dzieci nie było objętych ubezpieczeniem zdrowotnym, maksymalna liczba ubezpieczonych dzieci pozostających na utrzymaniu to 5, a najmniejsza 0.
smoker, liczba niepalących jest drastycznie większa niż liczba palących.
region, podział ze względu na pochodzenie (obszar) jest w miarę zrównoważony. Najwięcej pacjentów pochodzi z obszarów południowo-wschodnich.
charges, średni koszt ubezpieczenia medycznego kobiet wyniósł 12742.27 USD, a mężczyzn 14115.86 USD, największy wydatek był rzędu 63770.43 USD, a najmniejszy 1 135.94 USD. 50% wartości mieściło się w przedziale 5000.98 - 14489.39 USD (kobiety) i 4759.51 - 14437.11 USD. Odchylenie standardowe kosztów ubezpieczenia kobiet wynosi 11038.17 USD, a mężczyzn 13302.09 USD.

3.2 Wizualizacja zmiennych

Wizualizacja danych skupi się na identyfikacji wzorców kluczowych cech, które wyraźnie wyróżniają koszty leczenia pacjenta. W niektórych wizualizacjach posłużymy się podziałem ze względu na palaczy, w innych ze względu na płeć.

Na podstawie powyższego wykresu wydaje się, że wyższe opłaty są związane z palaczami w wieku starszym. Dodatkowo znacznie wyższe są odchylenia od linii regresji, co oznacza, że częsciej opłaty te są statystycznie wyższe w trakcie całego życia. Być może palacze bardziej zdają sobie sprawę z podatności na choroby. Może to być też zależność, że im starszy palacz, tym bardziej skomplikowane są jego warunki zdrowotne, które skutkują wyższymi opłatami za leczenie. Przy badaniu tych zależności zastosowaliśmy podział ze względu na palaczy, dla pewności zrobimy test ANOVA, który bada, czy istnieje różnica pomiędzy modelem bez interakcji z modelem zawierającym interakcję.

  Analysis of Variance Table
  
  Model 1: charges ~ age + smoker
  Model 2: charges ~ age * smoker
    Res.Df        RSS Df Sum of Sq      F  Pr(>F)  
  1    761 3.1590e+10                              
  2    760 3.1415e+10  1 174647951 4.2251 0.04017 *
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Na tym etapie można stwierdzić, że istnieje różnica pomiędzy modelem bez interakcji z modelem zawierającym interakcję, jednak wartość ta jest mało istotna statystycznie.

Na podstawie powyższego wykresu wydaje się, że wyższe opłaty są związane z palaczami z wyższym wskaźnikiem BMI. Być może im wyższe BMI palacza, tym bardziej jest on podatny na otyłość, co za tym idzie wyższe koszty leczenia.

  Analysis of Variance Table
  
  Model 1: charges ~ bmi + smoker
  Model 2: charges ~ bmi * smoker
    Res.Df        RSS Df  Sum of Sq     F    Pr(>F)    
  1    761 3.8518e+10                                  
  2    760 2.8169e+10  1 1.0349e+10 279.2 < 2.2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Można stwierdzić, że istnieje interakcja między zmienną smoker, a bmi.

Na podstawie powyższych wykresów można podejrzewać, że wyższe koszty leczenia są związane z mężczyznami. Może to być spowodowane tym, że mężczyźni są bardziej powszechnymi palaczami niż kobiety. Statystycznie uważa się także, że mężczyźni żyją krócej oraz mniej dbają o swoje zdrowie, co może skutkować wyższymi opłatami zdrowotnymi.

Sprawdzmy zależność pomiędzy opłatami, w poszczególnych regionach. Na podstawie powyższego wykresu można stwierdzić, że największe opłaty są ponoszone w regionie południowo wschodnim, co potwierdza wnioski analizy wzrokowej. Jednak warto sprawdzić, czy nie jest to wynikiem struktury osób palących oraz niepalących.

Na podstawie powyższego wykresu można stwierdzić istotną różnicę w rozkładzie palących w poszczególnych regionach, ze szczególnie ważną znaczną liczbą palaczy w regionie południowo-wschodnim. Może to oznaczać, że to nie region powoduje wyższe koszty ubezpieczeń. Przeprowadzona zostanie analiza ANOVA modeli region x opłaty z uwzględnieniem zmiennej palacz oraz jej braku. Dodatkowo przeprowadzona zostanie wizualizacja zależności opłaty x palacz:

Wizualizacja istotnie wskazuje na wyższe koszty leczenia związane z palaczami tytoniu. Osoby palące są bardziej podatne na choroby takie jak rak płuc, co skutkuje wyższymi opłatami za opiekę medyczną.

  Analysis of Variance Table
  
  Model 1: charges ~ region + smoker
  Model 2: charges ~ region * smoker
    Res.Df        RSS Df  Sum of Sq      F    Pr(>F)    
  1    759 4.2672e+10                                   
  2    756 4.0391e+10  3 2281107322 14.232 4.973e-09 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Powyższa analiza jednoznacznie potwierdza istotność modelu z uwzględnieniem palaczy, co oznacza, że region nie ma żadnego wpływu na zmiany w rozkładzie opłat. Można więc wyłączyć tę zmienną z dalszych analiz jako zmienną nieistotną. Natomiast można stwierdzić, że w regionie południowo wschodnim znajdzuje się najwięcej palaczy, co przekłada się na najwyższe opłaty w tym regionie.

Następnym krokiem jest przeanalizowanie liczby dzieci na wysokość opłat.

Jak można zauważyć na powyższym wykresie liczba dzieci nie ma wyraźnego związku z opłatami za ubezpieczenie. W związku z tym należy spodziewać się niskiej korelacji dla tej funkcji. Być może wysokie opłaty są niezależne od liczby dzieci. Dodatkowym zakłóceniem jest wysoki udział nieubezpieczonych dzieci, co może mieć wpływać na istotność tego czynnika danej populacji.

W tym przypadku będziemy przewidywać wartość zmiennej charges, czyli kosztów ubezpieczenia. Wykreślimy jej histogram.

Na podstawie powyższego histogramu możemy stwierdzić, że wysokość opłat koncentruje się przy minimalnych wartościach. Oznacza to, że większość ubezpieczeń opłacona jest najniższą stawką lub są one po prostu tanie.

3.3 Macierz korelacji

Za pomocą macierzy korelacji możemy wizualizować związek między zmienną zależną a zmiennymi niezależnymi, a także między nimi. Należy podkreślić, że dotyczy to tylko zmiennych liczbowych. Potwierdzi to więc poszczególne wnioski dla liczbowych zależności, które zostały uprzednio podjęte:

Te same wartości korelacji możemy zobrazować przy pomocy mapy kolorów, gdzie im cieplejsza barwa, tym silniejsza korelacja między zmiennymi numerycznymi.

4 Wstępne spostrzeżenia

Eksploracja danych ukierunkowała na następujące wnioski:

Do wyższych kosztów ubezpieczenia prowadzą:
- starszy wiek,
- palenie tytoniu,
- płeć - mężczyzna,
- bmi na poziomie > 30.
Na podstawie segmentu korelacji macierzy, age i bmi ma największą korelację z charges
Istnieje interakcja między smoker, a bmi.

5 Propozycje modeli nauczone na podstawie danych uczących

Przedstawimy zestaw różnych propozycji modeli nauczonych na podstawie danych uczących. Podsumujemy każdy z modeli wyjaśniając wyniki wyjściowe działania funkcji summary;

Sekcja Residuals zawiera podstawowe informacje o resztach modelu (mediana, minimum, maksimum).
Sekcja Coefficients zawiera znalezione współczynniki modelu ich błędy standardowe, informacje o ich istotności statystycznej oraz wartości statystyki \(t\).
Residual standard error, znajdujący się w ostatniej sekcji wyników regresji, jest estymatorem parametru \(\sigma\) (czyli wariancji błędu).
F-statistic znajdująca się u dołu tabeli testuje hipotezę zerową, iż nasz model składa się wyłącznie ze stałego czynnika (jest to tak zwany model zerowy), tzn. zmienna objaśniana nie zależy od zmiennej objaśniającej.
Kwadrat współczynnika korelacji \(r^2\) nazywany jest często współczynnikiem determinacji. Określa on jaka część całkowitej sumy kwadratów jest wyjaśniana przez model.
U dołu tabeli znajduje się również zmodyfikowany współczynnik determinacji (Adjusted R-squared). Jak możemy zauważyć przyjmuje on wartość (niemalże) taką samą jak \(r^2\).

Wyliczymy błąd dla każdego ze stworzonych modeli oraz zobrazujemy wykresy diagnostyczne, przy pomocy funkcji plot;

Pierwszy wykres przedstawia reszty wykreślone względem wartości dopasowanych. Punkty powinny być równomiernie rozrzucone względem osi x. Czerwona linia (wygładzenie) wskazuje na lekkie niedopasowane modelu dla mniejszych wartości.
Drugim wykresem jest wykres kwantyl-kwantyl standaryzowanych reszt. Reszty standaryzowane, są to reszty podzielone przez odchylenie standardowe. Wykres kwantyl-kwantyl powinien przypominać linię prostą. Pewne niedopasowanie można zauważyć w przypadku bardziej ekstremalnych obserwacji.
Trzeci wykres jest wygodny do badania homoskedastyczności (czerwona krzywa powinna być jak najbardziej pozioma).
Czwarty wykres służy do badania, które obserwacje są wpływowe. Uwidoczniona jest na nim odległość Cooka, która bada różnicę predykcji modelu wyjściowego z modelem po usunięciu obserwacji.

Model 1 Zaczniemy od najprostszego modelu, gdzie predyktorem będzie najbardziej skorelowana zmienna, czyli age.

  
  Call:
  lm(formula = charges ~ age, data = insurance_train)
  
  Residuals:
     Min     1Q Median     3Q    Max 
   -8054  -6582  -5907   5120  47815 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)   3007.5     1252.8   2.401   0.0166 *  
  age            261.6       29.6   8.838   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 11660 on 762 degrees of freedom
  Multiple R-squared:  0.09297, Adjusted R-squared:  0.09178 
  F-statistic:  78.1 on 1 and 762 DF,  p-value: < 2.2e-16

Błąd modelu:

  [1] 103681190558

Wykresy diagnostyczne:

Model 2 Zobaczmy, jak wygląda nasz model, gdy zmienna zależna zostanie zlogarytmowana.

  
  Call:
  lm(formula = log(charges) ~ age, data = insurance_train)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -1.3455 -0.4006 -0.2997  0.4382  2.1873 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 7.771840   0.082592   94.10   <2e-16 ***
  age         0.033826   0.001951   17.34   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.769 on 762 degrees of freedom
  Multiple R-squared:  0.2829,  Adjusted R-squared:  0.2819 
  F-statistic: 300.6 on 1 and 762 DF,  p-value: < 2.2e-16

W porównaniu do wcześniejszego modelu, widzimy poprawę wartości Adjusted R-squared oraz spadek błędu standardowego reszt.

Błąd modelu:

  [1] 450.5896

Wykresy diagnostyczne:

Model 3 Rozważmy model \(\log(\mathrm{charges}) = a + b \cdot \mathrm{age} + c \cdot \mathrm{age}^2\).

  
  Call:
  lm(formula = log(charges) ~ age + I(age^2), data = insurance_train)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -1.2410 -0.4451 -0.2888  0.4950  2.2897 
  
  Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
  (Intercept)  7.2444242  0.2356179  30.746  < 2e-16 ***
  age          0.0641783  0.0128526   4.993 7.36e-07 ***
  I(age^2)    -0.0003808  0.0001594  -2.389   0.0171 *  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.7666 on 761 degrees of freedom
  Multiple R-squared:  0.2882,  Adjusted R-squared:  0.2863 
  F-statistic: 154.1 on 2 and 761 DF,  p-value: < 2.2e-16

Brak poprawy wartości Adjusted R-squared, ta sama tendencja jest również zauważalna przy błędzie standardowym reszt.

Błąd modelu:

  [1] 447.2353

Wykresy diagnostyczne:

Model 4 Zazwyczaj staramy się wyjaśnić wartości zmiennej zależnej za pomocą więcej niż jednego predyktora. Mówimy wtedy, że mamy do czynienia z regresją wieloraką. Drugą pod względem wielkości korelację z charges miało bmi, dodajmy więc tą zmienną do naszego pierwotnego modelu.

  
  Call:
  lm(formula = charges ~ age + bmi, data = insurance_train)
  
  Residuals:
     Min     1Q Median     3Q    Max 
  -14394  -7001  -5064   4864  48078 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept) -7060.22    2339.86  -3.017  0.00263 ** 
  age           249.19      29.23   8.525  < 2e-16 ***
  bmi           341.27      67.41   5.063 5.19e-07 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 11480 on 761 degrees of freedom
  Multiple R-squared:  0.1225,  Adjusted R-squared:  0.1202 
  F-statistic: 53.13 on 2 and 761 DF,  p-value: < 2.2e-16

Dostrzegamy małą wartość Adjusted R-squared, dodatkowo błąd standardowy reszt jest wysoki.

Błąd modelu:

  [1] 100302825875

Wykresy diagnostyczne:

Model 5 Zmodyfikujmy model nr 4 poprzez zlogarytmowanie zmiennej zależnej.

  
  Call:
  lm(formula = log(charges) ~ age + bmi, data = insurance_train)
  
  Residuals:
      Min      1Q  Median      3Q     Max 
  -1.5694 -0.4325 -0.2800  0.4689  2.0766 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept) 7.453635   0.156239  47.706   <2e-16 ***
  age         0.033435   0.001952  17.130   <2e-16 ***
  bmi         0.010786   0.004501   2.396   0.0168 *  
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 0.7666 on 761 degrees of freedom
  Multiple R-squared:  0.2882,  Adjusted R-squared:  0.2864 
  F-statistic: 154.1 on 2 and 761 DF,  p-value: < 2.2e-16

Poprawa wartości Adjusted R-squared oraz spadek błędu standardowego reszt.

Błąd modelu:

  [1] 447.2148

Wykresy diagnostyczne:

Model 6 Stworzymy model, gdzie zmienną numeryczną przewiduję za pomocą zmiennej kategorycznej. Mówimy wtedy, że mamy do czynienia z analizą wariancji (ANOVA).

  
  Call:
  lm(formula = charges ~ smoker, data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -19685.4  -4892.7   -575.7   3699.1  31255.6 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept)   8536.6      304.4   28.04   <2e-16 ***
  smokeryes    23978.2      673.6   35.59   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 7506 on 762 degrees of freedom
  Multiple R-squared:  0.6244,  Adjusted R-squared:  0.6239 
  F-statistic:  1267 on 1 and 762 DF,  p-value: < 2.2e-16

Widzimy ogromną poprawę Adjusted R-squared, w porównaniu do wszystkich poprzednich modeli.

Błąd modelu:

  [1] 42929636503

Wykresy diagnostyczne:

Model 7 Zmienną numeryczną możemy także wyjaśnić za pomocą innej zmiennej numerycznej i zmiennej kategorycznej jednocześnie. Jest to sytuacja która jest połączeniem regresji wielorakiej i analizy wariancji. W takim przypadku przyjęło się używać nazwy analiza kowariancji (ANCOVA).

  
  Call:
  lm(formula = charges ~ age + smoker, data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -16288.8  -1999.1  -1265.9   -229.7  28591.9 
  
  Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
  (Intercept) -2258.01     703.44   -3.21  0.00138 ** 
  age           270.21      16.35   16.53  < 2e-16 ***
  smokeryes   24099.55     578.30   41.67  < 2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 6443 on 761 degrees of freedom
  Multiple R-squared:  0.7236,  Adjusted R-squared:  0.7229 
  F-statistic: 996.3 on 2 and 761 DF,  p-value: < 2.2e-16

Błąd modelu:

  [1] 31589934284

Wykresy diagnostyczne:

Model 8 Dopasujmy teraz model, który uwzględnia również trzecią zmienną, mianowicie bmi.

  
  Call:
  lm(formula = charges ~ age + bmi + smoker, data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -12358.8  -3121.0   -955.4   1530.5  28859.7 
  
  Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
  (Intercept) -12254.10    1248.36  -9.816   <2e-16 ***
  age            257.92      15.53  16.611   <2e-16 ***
  bmi            338.90      35.80   9.466   <2e-16 ***
  smokeryes    24091.75     547.31  44.018   <2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 6098 on 760 degrees of freedom
  Multiple R-squared:  0.7528,  Adjusted R-squared:  0.7518 
  F-statistic: 771.4 on 3 and 760 DF,  p-value: < 2.2e-16

Trzecia zmienna, bmi, spowodowała dość istotną zmianę wartości błędu standardowego reszt. Podobną zależność widzimy w wartości współczynnika determinacji, czyli nasz nowy model wyjaśnia większy procent wariancji kosztu leczenia/ubezpieczenia.

Błąd modelu:

  [1] 28258366146

Wykresy diagnostyczne:

Model 9 Interakcje pozwalają lepiej modelować zależności pomiędzy predyktorami. Wiemy, że koszt ubezpieczenia wzrasta wraz ze wskaźnikiem bmi. Wiemy, też że wysokość kosztów leczenia zależy od tego czy pacjent jest palaczem, czy nie. Odniesiemy się do tej kwestii, dopasowując model z interakcjami.

  
  Call:
  lm(formula = charges ~ bmi * smoker, data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -19827.7  -4388.0   -614.5   2878.9  30946.9 
  
  Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
  (Intercept)     6640.53    1281.95   5.180 2.84e-07 ***
  bmi               61.27      40.65   1.507    0.132    
  smokeryes     -19650.61    2667.44  -7.367 4.56e-13 ***
  bmi:smokeryes   1409.61      84.36  16.709  < 2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 6088 on 760 degrees of freedom
  Multiple R-squared:  0.7536,  Adjusted R-squared:  0.7526 
  F-statistic: 774.7 on 3 and 760 DF,  p-value: < 2.2e-16

Błąd modelu:

  [1] 28169423596

Wykresy diagnostyczne:

Model 10 Warto także zbudować model, który włącza wszystkie dostępne zmienne.

  
  Call:
  lm(formula = charges ~ ., data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -11441.6  -2894.1   -898.4   1413.3  29532.5 
  
  Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
  (Intercept)     -12540.27    1302.85  -9.625  < 2e-16 ***
  age                257.41      15.50  16.605  < 2e-16 ***
  sexmale            -72.81     441.99  -0.165  0.86919    
  bmi                338.05      37.77   8.950  < 2e-16 ***
  children           506.61     179.07   2.829  0.00479 ** 
  smokeryes        24087.00     548.81  43.890  < 2e-16 ***
  regionnorthwest    -21.07     630.01  -0.033  0.97333    
  regionsoutheast   -259.91     641.74  -0.405  0.68558    
  regionsouthwest   -517.30     638.75  -0.810  0.41828    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 6082 on 755 degrees of freedom
  Multiple R-squared:  0.7557,  Adjusted R-squared:  0.7531 
  F-statistic: 291.9 on 8 and 755 DF,  p-value: < 2.2e-16

Jak możemy zauważyć, nie wszystkie zastosowane zmienne są istotne statystycznie.

Błąd modelu:

  [1] 27928236437

Wykresy diagnostyczne:

Model 11 Wiemy, że wraz ze wzrostem wieku koszty leczenia są większe. Aby uchwycić tę nieliniową zależność, do modelu zawierającego wszystkie zmienne dodamy nową zmienną \(age^{2}\).

  
  Call:
  lm(formula = charges ~ age + I(age^2) + children + bmi + sex + 
      smoker + region, data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -11511.1  -2928.5   -884.5   1333.9  30512.7 
  
  Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
  (Intercept)     -6524.963   2213.723  -2.948 0.003302 ** 
  age               -94.291    106.133  -0.888 0.374595    
  I(age^2)            4.409      1.316   3.349 0.000851 ***
  children          691.952    186.284   3.714 0.000219 ***
  bmi               334.139     37.534   8.902  < 2e-16 ***
  sexmale          -113.944    439.198  -0.259 0.795370    
  smokeryes       24087.048    545.130  44.186  < 2e-16 ***
  regionnorthwest    -5.206    625.811  -0.008 0.993365    
  regionsoutheast  -197.266    637.720  -0.309 0.757156    
  regionsouthwest  -457.334    634.726  -0.721 0.471427    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 6041 on 754 degrees of freedom
  Multiple R-squared:  0.7593,  Adjusted R-squared:  0.7564 
  F-statistic: 264.2 on 9 and 754 DF,  p-value: < 2.2e-16

Błąd modelu:

  [1] 27518849319

Wykresy diagnostyczn:

Model 12 Zmienimy zmienną BMI na binarną, aby wskazać stan otyłości (BMI >= 30). Chcemy uchwycić związek między otyłością a wyższymi kosztami leczenia.

  
  Call:
  lm(formula = charges ~ age + I(age^2) + children + bmi + bmi30 + 
      sex + smoker + region, data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -12346.3  -3307.0     89.2   1507.2  29063.9 
  
  Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
  (Intercept)     -3121.865   2401.538  -1.300 0.194018    
  age               -75.037    105.487  -0.711 0.477094    
  I(age^2)            4.164      1.308   3.183 0.001520 ** 
  children          660.691    185.114   3.569 0.000381 ***
  bmi               172.113     59.307   2.902 0.003815 ** 
  bmi30            2514.267    716.064   3.511 0.000473 ***
  sexmale          -122.815    435.943  -0.282 0.778234    
  smokeryes       24124.604    541.186  44.577  < 2e-16 ***
  regionnorthwest  -156.107    622.647  -0.251 0.802103    
  regionsoutheast  -146.791    633.146  -0.232 0.816722    
  regionsouthwest  -551.643    630.583  -0.875 0.381954    
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 5996 on 753 degrees of freedom
  Multiple R-squared:  0.7631,  Adjusted R-squared:   0.76 
  F-statistic: 242.6 on 10 and 753 DF,  p-value: < 2.2e-16

Błąd modelu:

  [1] 27075545022

Wykresy diagnostyczne:

Model 13 Załóżmy, że sensowne jest przetestowanie połączonego efektu dwóch zmiennych. Zbudujemy model z interakcją między zmienną bmi30 a smoker.

  
  Call:
  lm(formula = charges ~ age + I(age^2) + children + bmi + bmi30 * 
      smoker + sex + region, data = insurance_train)
  
  Residuals:
       Min       1Q   Median       3Q      Max 
  -17446.0  -1605.2  -1123.4   -632.3  24014.6 
  
  Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
  (Intercept)       686.5582  1771.8680   0.387  0.69851    
  age               -81.2339    77.5487  -1.048  0.29520    
  I(age^2)            4.2588     0.9619   4.427 1.10e-05 ***
  children          737.1084   136.1190   5.415 8.24e-08 ***
  bmi               129.0756    43.6322   2.958  0.00319 ** 
  bmi30           -1115.1845   545.5704  -2.044  0.04129 *  
  smokeryes       13509.4030   577.9168  23.376  < 2e-16 ***
  sexmale          -633.6407   321.1148  -1.973  0.04883 *  
  regionnorthwest  -290.8889   457.7658  -0.635  0.52533    
  regionsoutheast  -601.6125   465.7994  -1.292  0.19690    
  regionsouthwest  -996.6875   463.9016  -2.148  0.03199 *  
  bmi30:smokeryes 20166.5375   796.3313  25.324  < 2e-16 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  Residual standard error: 4408 on 752 degrees of freedom
  Multiple R-squared:  0.8722,  Adjusted R-squared:  0.8703 
  F-statistic: 466.4 on 11 and 752 DF,  p-value: < 2.2e-16

Jak dotychczas, w przypadku tego modelu, uzyskaliśmy największą wartość Adjusted R-squared.

Błąd modelu:

  [1] 14613156386

Wykresy diagnostyczne:

6 Wybór najlepszego modelu predykcyjnego

Jak pokazaliśmy wyżej, możemy dopasować wiele modeli wyjaśniających zmienność kosztów leczenia. W zależności od zastosowań zasadne jest używanie różnych kryteriów wyboru modelu. Najlepsze modele wybierzemy na podstawie weryfikacji na danych testowych i ogólnej wiedzy statystycznej oraz odpowiednich testach.

Model nr 1

Tworzymy nową zmienną - wartość przewidywaną dla zmiennej charges.

insurance_test$pred <- predict(model_1, insurance_test)

Sprawdzamy, jak dobrze przewidywane wartości są skorelowane z tymi rzeczywistymi.

  [1] 0.2897713

Wyliczyliśmy ważny współczynnik oceny prognozy, tj. błąd średni kwadratowy (MSE).

  [1] 130429983

Model nr 2