EDA: Badanie czynników kształtujących zdrowie

Wprowadzenie

Zdrowie to znacznie więcej niż brak chorób – stanowi fundament codziennego dobrostanu oraz poziomu energii, z jakim funkcjonujemy w społeczeństwie. Jak trafnie zauważył Denis Waitley:
„Czas i zdrowie to dwa cenne skarby, których nie rozpoznajemy i nie doceniamy, dopóki się nie wyczerpią.”

Dzięki analizie danych dotyczących diety i ruchu, podzielonych na kategorie demograficzne, sprawdzimy, w jakim stopniu czynniki zewnętrzne takie jak status materialny czy wykształcenie – modelują nasze podejście do dbałości o siebie.

Zmienne

State - nazwa stanu
Class - główna kategoria tematyczna wskaźnika (np. Obesity / Nutrition)
Question - pytanie zadane w badaniu
Data_Value -oszacowany odsetek (%)
Low_Confidence_Limit - dolna granica przedziału ufności
High_Confidence_Limit - górna granica przedziału ufności
Sample_Size - rozmiar badanej próbki
Strat_Category - Kategoria (np. “Income”, “Education,”Race/Ethnicity”)
Strat_Group - podgrupa w kategorii (np. “18-24”, “Male”, “15,000 - 24,999”).

Pakiety, których będę używał

library(readxl)
library(dplyr)
library(tidyr)
library(mice)
library(ggplot2)
library(corrplot)

Ogląd oraz czyszczenie danych

Dane pochodzą ze strony Data.Gov: Nutrition, Physical Activity, and Obesity - Behavioral Risk Factor Surveillance System [https://catalog.data.gov/dataset/nutrition-physical-activity-and-obesity-behavioral-risk-factor-surveillance-system]

dane <- read_excel("data/Nutrition_PhysicalActivity_Obesity.xlsx", sheet = "Analyzed_Data")

str(dane)

## tibble [110,880 × 10] (S3: tbl_df/tbl/data.frame)
##  $ YearStart              : num [1:110880] 2011 2011 2011 2011 2011 ...
##  $ LocationDesc           : chr [1:110880] "Alabama" "Alabama" "Alabama" "Alabama" ...
##  $ Class                  : chr [1:110880] "Obesity / Weight Status" "Obesity / Weight Status" "Obesity / Weight Status" "Obesity / Weight Status" ...
##  $ Question               : chr [1:110880] "Percent of adults aged 18 years and older who have obesity" "Percent of adults aged 18 years and older who have obesity" "Percent of adults aged 18 years and older who have obesity" "Percent of adults aged 18 years and older who have obesity" ...
##  $ Data_Value             : chr [1:110880] "34.8" "35.8" "32.3" "34.1" ...
##  $ Low_Confidence_Limit   : chr [1:110880] "31.3" "31.1" "28.0" "29.7" ...
##  $ High_Confidence_Limit  : chr [1:110880] "38.5" "40.8" "36.8" "38.8" ...
##  $ Sample_Size            : num [1:110880] 1367 757 861 785 1125 ...
##  $ StratificationCategory1: chr [1:110880] "Income" "Income" "Income" "Income" ...
##  $ Stratification1        : chr [1:110880] "$15,000 - $24,999" "$25,000 - $34,999" "$35,000 - $49,999" "$50,000 - $74,999" ...

summary(dane)

##    YearStart    LocationDesc          Class             Question        
##  Min.   :2011   Length:110880      Length:110880      Length:110880     
##  1st Qu.:2014   Class :character   Class :character   Class :character  
##  Median :2017   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :2017                                                           
##  3rd Qu.:2020                                                           
##  Max.   :2024                                                           
##                                                                         
##   Data_Value        Low_Confidence_Limit High_Confidence_Limit  Sample_Size    
##  Length:110880      Length:110880        Length:110880         Min.   :    50  
##  Class :character   Class :character     Class :character      1st Qu.:   494  
##  Mode  :character   Mode  :character     Mode  :character      Median :  1079  
##                                                                Mean   :  3625  
##                                                                3rd Qu.:  2399  
##                                                                Max.   :476876  
##                                                                NA's   :13214   
##  StratificationCategory1 Stratification1   
##  Length:110880           Length:110880     
##  Class :character        Class :character  
##  Mode  :character        Mode  :character  
##                                            
##                                            
##                                            
##

Dla czytelności zmieniam nazwy kolumn

dane <- dane %>% 
  rename(Year = YearStart, 
         State = LocationDesc ,
         Strat_Category = StratificationCategory1,
         Strat_Group = Stratification1
         )

Widzę, że pozycje Data_Value, Low_Confidence_Limit oraz High_Confidence_Limit nie są numeryczne.
Zmieniam:

dane <- dane %>%
  mutate(
    Data_Value = as.numeric(Data_Value),
    Low_Confidence_Limit = as.numeric(Low_Confidence_Limit),
    High_Confidence_Limit = as.numeric(High_Confidence_Limit)
  )

str(dane)

## tibble [110,880 × 10] (S3: tbl_df/tbl/data.frame)
##  $ Year                 : num [1:110880] 2011 2011 2011 2011 2011 ...
##  $ State                : chr [1:110880] "Alabama" "Alabama" "Alabama" "Alabama" ...
##  $ Class                : chr [1:110880] "Obesity / Weight Status" "Obesity / Weight Status" "Obesity / Weight Status" "Obesity / Weight Status" ...
##  $ Question             : chr [1:110880] "Percent of adults aged 18 years and older who have obesity" "Percent of adults aged 18 years and older who have obesity" "Percent of adults aged 18 years and older who have obesity" "Percent of adults aged 18 years and older who have obesity" ...
##  $ Data_Value           : num [1:110880] 34.8 35.8 32.3 34.1 28.8 16.3 27.8 35.2 35.5 38 ...
##  $ Low_Confidence_Limit : num [1:110880] 31.3 31.1 28 29.7 25.4 12.6 14.4 30.7 31.6 34.5 ...
##  $ High_Confidence_Limit: num [1:110880] 38.5 40.8 36.8 38.8 32.5 20.9 46.9 40 39.6 41.5 ...
##  $ Sample_Size          : num [1:110880] 1367 757 861 785 1125 ...
##  $ Strat_Category       : chr [1:110880] "Income" "Income" "Income" "Income" ...
##  $ Strat_Group          : chr [1:110880] "$15,000 - $24,999" "$25,000 - $34,999" "$35,000 - $49,999" "$50,000 - $74,999" ...

Do dalszej analizy warto będzie przypisać stany do 4 głównych regionów USA (według klasyfikacji Census)
Można zobaczyć jak to wygląda na mapie: [https://www2.census.gov/geo/pdfs/maps-data/maps/reference/us_regdiv.pdf]

state_map <- data.frame(
  State = state.name,
  Region = state.region
  )

Upewnijmy się, że wszystkie dane są ze sobą zgodne

zgranie <- setdiff(unique(dane$State), state_map$State)

zgranie

## [1] "District of Columbia" "National"             "Guam"                
## [4] "Puerto Rico"          "Virgin Islands"

Pojawiła mi się nieoczekiwana zmienna “National” (oprócz 4 terytoriów, które nie są uznawane za stany i nie znajdują się w state_map)

dane %>% 
  filter(State == "National") %>% 
  summary()

##       Year         State              Class             Question        
##  Min.   :2011   Length:2016        Length:2016        Length:2016       
##  1st Qu.:2014   Class :character   Class :character   Class :character  
##  Median :2017   Mode  :character   Mode  :character   Mode  :character  
##  Mean   :2017                                                           
##  3rd Qu.:2020                                                           
##  Max.   :2024                                                           
##                                                                         
##    Data_Value    Low_Confidence_Limit High_Confidence_Limit  Sample_Size    
##  Min.   : 8.70   Min.   : 7.50        Min.   :10.10         Min.   :   664  
##  1st Qu.:25.80   1st Qu.:24.70        1st Qu.:26.90         1st Qu.: 26336  
##  Median :32.00   Median :30.80        Median :33.10         Median : 50080  
##  Mean   :31.83   Mean   :30.69        Mean   :33.02         Mean   : 87382  
##  3rd Qu.:36.60   3rd Qu.:35.50        3rd Qu.:37.80         3rd Qu.:117780  
##  Max.   :69.30   Max.   :68.80        Max.   :69.80         Max.   :476876  
##  NA's   :3       NA's   :3            NA's   :3             NA's   :3       
##  Strat_Category     Strat_Group       
##  Length:2016        Length:2016       
##  Class :character   Class :character  
##  Mode  :character   Mode  :character  
##                                       
##                                       
##                                       
##

Odkryliśmy tutaj przypadkiem bardzo ważną rzecz - w naszych danych istnieje wynik zbiorczy dla całych Stanów Zjednoczonych (ogólnokrajowa estymacja). Gdybym to zostawił to w analizach porównujących stany byłyby dodatkową obserwacją i zniekształcałyby statystyki stanowe. Usuńmy tę pozycje z głównego zestawu (zostawimy ją w innej zmiennej by ewentualnie móc porównywać wyniki w stanach do poziomu krajowego)

dane_national <- dane %>%
  filter(State == "National")

dane <- dane %>%
  filter(State != "National")

Z analizy wykluczę obszary nieprzypisane do żadnego regionu (poza DC), ponieważ pochodzące z nich dane mogą być niewiarygodne. Wyjątek stanowi District of Columbia – ze względu na status stolicy administracyjnej oraz bezpośrednie sąsiedztwo z regionem South, zostanie on włączony do badanego zbioru.

dane_state <- dane %>% 
  filter(!(State %in% c("Guam", "Puerto Rico", "Virgin Islands")))

state_map <- state_map %>%  
  add_row(State = "District of Columbia", Region = "South")

I teraz połączmy tabelę state_map (czyli stan + region) z naszym danymi z samymi stanami

dane_state <- dane_state %>%
  left_join(state_map, by = "State")

Dla niektórych zmiennych lepiej działać na typie factor

dane_state <- dane_state %>%
  mutate(
    State = factor(State),
    Region = factor(Region),
    Class = factor(Class),
    Strat_Category = factor(Strat_Category),
    Strat_Group = factor(Strat_Group)
  )

Co więcej niektóre grupy powinny być uszeregowane (np dochód czy wiek). Jednak z uwagi, że kolumna Strat_Group posiada wszystkie kategorie to problematyczne byłoby nadawanie poziomów dla ich wszystkich na raz. Zajmiemy się tym, gdy będzie to potrzebne w konkretnej analizie.

Sprawdźmy jak wygląda sytuacja z wartościami pustymi

sum(is.na(dane_state))

## [1] 40904

Braki występują w blisko 40% danych (40 904/102 816). Aby przyjrzeć się temu problemowi bliżej, wykorzystam pakiet mice

md.pattern(dane_state, rotate.names = TRUE)

##       Year State Class Question Strat_Category Strat_Group Region Data_Value
## 92590    1     1     1        1              1           1      1          1
## 10226    1     1     1        1              1           1      1          0
##          0     0     0        0              0           0      0      10226
##       Low_Confidence_Limit High_Confidence_Limit Sample_Size      
## 92590                    1                     1           1     0
## 10226                    0                     0           0     4
##                      10226                 10226       10226 40904

Wszystko jasne - skoro mamy braki w kolumnie Data_Value to również będą występować braki w kolumnach (LCL,HCL, SC, SG)
Zobaczmy jak to jest z kategoriami

braki_kategorie <- dane_state %>% 
  filter(is.na(Data_Value)) %>% 
  count(Strat_Category, sort = TRUE)

knitr::kable(braki_kategorie, caption = "Braki danych w poszczególnych kategoriach")

Braki danych w poszczególnych kategoriach
Strat_Category	n
Race/Ethnicity	9601
Income	222
Age (years)	186
Education	124
Sex	62
Total	31

Największa liczba braków danych występuje w kategorii Race/Ethnicity. Może to wynikać z faktu, że w wielu stanach liczebność konkretnych mniejszości etnicznych (np. nielicznej ilości Azjatów mieszkających w Alabamie) jest zbyt niska, aby uzyskać reprezentatywną próbę. W takich przypadkach dane są często celowo pomijane w celu uniknięcia wyciągania wniosków na podstawie zbyt małej grupy badawczej. Jednak na ten moment nie będę usuwał rekordów związanych z tą kategorią.

Chciałbym się upewnić, czy na przestrzeni lat badacze konsekwentnie zbierali te same dane. Być może zauważono, że w niektórych stanach (jak we wspomnianej Alabamie) liczebność pewnych grup jest stale zbyt niska, co mogło wpłynąć na decyzję o zmianie sposobu raportowania w kolejnych latach i spowodować braki w danych.

c <- table(dane_state$Strat_Group)
knitr::kable(c, col.names = c("Grupa", "Liczebność"), caption = "Podsumowanie grup")

Podsumowanie grup
Grupa	Liczebność
$15,000 - $24,999	3672
$25,000 - $34,999	3672
$35,000 - $49,999	3672
$50,000 - $74,999	3672
$75,000 or greater	3672
18 - 24	3672
2 or more races	3672
25 - 34	3672
35 - 44	3672
45 - 54	3672
55 - 64	3672
65 or older	3672
American Indian/Alaska Native	3672
Asian	3672
College graduate	3672
Data not reported	3672
Female	3672
Hawaiian/Pacific Islander	3672
High school graduate	3672
Hispanic	3672
Less than $15,000	3672
Less than high school	3672
Male	3672
Non-Hispanic Black	3672
Non-Hispanic White	3672
Other	3672
Some college or technical school	3672
Total	3672

Dane są bardzo dobrze zebrane – każda grupa stratyfikacji występuje w danych dokładnie tyle samo razy (3672 obserwacje na grupę). Takie coś pokazuje, że badacze przyjęli stały schemat raportowania, w którym każdy stan dostarcza informacje o wszystkich grupach w stałych odstępach czasu.

Żeby nie było nieporozumień: ewentualne braki danych (NA), o których wspomniałem wcześniej, nie wynikają z pominięcia wierszy, lecz z braku konkretnych wartości dla już istniejących rekordów.

Liczba 3672 zapewne wzięła się z prostego iloczynu: $lata * ilość stanów * ilość pytań$
Upewnijmy się:

iloczyn<- length(unique(dane_state$Year)) * length(unique(dane_state$State)) * length(unique(dane_state$Question))

iloczyn

## [1] 6426

Iloczyn wyszedł większy niż wcześniej ukazana wartość 3672 (brakuje nam 2754). Skoro dla każdego roku mamy 50 stanów (w tym DC, wiec ogólnie 51), to być może nie zawsze zadawane było każde pytanie…

Zobaczmy ile razy każde pytanie występuje w każdym roku

mapa_pytan <- table(dane_state$Question, dane_state$Year)

knitr::kable(mapa_pytan, caption = "Pytania zadawane na przestrzeni lat")

Pytania zadawane na przestrzeni lat
	2011	2012	2013	2014	2015	2016	2017	2018	2019	2020	2021	2022	2023	2024
Percent of adults aged 18 years and older who have an overweight classification	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428
Percent of adults aged 18 years and older who have obesity	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	1428	0	1428	0	1428	0	1428	0	1428	0	0	0	1428	0
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic physical activity (or an equivalent combination) and engage in muscle-strengthening activities on 2 or more days a week	1428	0	1428	0	1428	0	1428	0	1428	0	0	0	1428	0
Percent of adults who achieve more than 300 minutes a week of moderate-intensity aerobic physical activity or 150 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	1428	0	1428	0	1428	0	1428	0	1428	0	0	0	1428	0
Percent of adults who engage in muscle-strengthening activities on 2 or more days a week	1428	0	1428	0	1428	0	1428	0	1428	0	0	0	1428	0
Percent of adults who engage in no leisure-time physical activity	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428	1428
Percent of adults who report consuming fruit less than one time daily	0	0	0	0	0	0	1428	0	1428	0	1428	0	0	0
Percent of adults who report consuming vegetables less than one time daily	0	0	0	0	0	0	1428	0	1428	0	1428	0	0	0

Moje wątpliwości zostały wyjaśnione. Nie w każdym roku grupa badawcza zadawała uczestnikom te same pytania.
Jednocześnie liczba 1428 się zgadza: $51 jednostek * 28 grup$.
Czyli jeżeli w danym roku zdecydowano się zadać dane pytanie to zostało ono zebrane dla wszystkich 51 jednostek z podziałem na nasze 28 grup.

Ile dokładnie mamy braków?

table(mapa_pytan)

## mapa_pytan
##    0 1428 
##   54   72

I tu się rozwiązuje problem z brakującymi 2754 kombinacjami.

Faktycznie zadanych pytań przez 14 lat: $72$
Oraz 51 stanów daje: $51 * 72 = 3672$
Jak i nasze braki: $54 * 51 = 2754$

Dla lepszego rozeznania zobaczmy jak to wygląda na wykresie:
(Z uwagi na długość pytań, które mogą wpłynąć na czytelność wykresu przypisze im konkretne ID)

legenda_pytan <- dane_state %>%
  distinct(Question) %>%
  mutate(ID = paste0("P", row_number()))

mapa_danych <- dane_state %>%
  distinct(Year, Question) %>%
  left_join(legenda_pytan, by = "Question")

mapa_danych %>% 
  ggplot(aes(x = as.factor(Year), y = ID)) +
  geom_tile(fill = "steelblue3", color = "white", linewidth = 0.8) + 
  labs(
    title = "Dostępność danych w czasie",
    x = "Rok",
    y = "Numer pytania"
  ) +
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

knitr::kable(legenda_pytan, caption = "Legenda Pytań")

Legenda Pytań
Question	ID
Percent of adults aged 18 years and older who have obesity	P1
Percent of adults aged 18 years and older who have an overweight classification	P2
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	P3
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic physical activity (or an equivalent combination) and engage in muscle-strengthening activities on 2 or more days a week	P4
Percent of adults who achieve more than 300 minutes a week of moderate-intensity aerobic physical activity or 150 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	P5
Percent of adults who engage in muscle-strengthening activities on 2 or more days a week	P6
Percent of adults who engage in no leisure-time physical activity	P7
Percent of adults who report consuming fruit less than one time daily	P8
Percent of adults who report consuming vegetables less than one time daily	P9

Teraz już widzimy, jak wygląda kwestia z pytaniami w poszczególnych latach - jedynie 3 pytania były sukcesywnie zadawane z roku na rok. Tym samym rozwiązała się kwestia braku danych. Moim zdaniem nie trzeba usuwać tych pustych miejsc - będziemy analizować zależności w różnych przedziałach czasowych.

Przejdźmy do kolejnych zmiennych - granic przedziału ufności.
Wpierw zobaczmy czy istnieje sytuacja gdzie Data_Value nie mieści się w przedziale ufności

dane_state %>% 
  filter(Data_Value < Low_Confidence_Limit | Data_Value > High_Confidence_Limit) %>% 
  nrow()

## [1] 0

Bardzo dobra informacja. We wszystkich rekordach wartość estymowana (Data_Value) mieści się w raportowanych przedziałach ufności. Gdybyśmy otrzymali wartość różną od zera to trzeba by było sprawdzić czy dany rekord nie zawiera błędu (przedział ufności jest budowany wokół wartości estymowanej)

Popatrzmy na podstawowe statystyki związane z szerokością przedziału ufności:

summary(dane_state$High_Confidence_Limit -dane_state$Low_Confidence_Limit )

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   1.200   5.100   7.500   9.369  11.000  59.200   10226

Połowa naszych danych (czyli ponad 55 tys.) ma widełki nie szersze niż 7.5 punktu procentowego - to bardzo dobra wiadomość.
Uwagę od razu przykuwa wartość maksymalna 59.2 - jest to bardzo duży przedział ufności i tak naprawdę wynik tego badania może być mało precyzyjne. Zidentyfikujmy ten rekord:

dane_state %>% 
  filter(High_Confidence_Limit -Low_Confidence_Limit > 59) %>% 
  select(Year, State, Class, Strat_Category, Data_Value)

## # A tibble: 1 × 5
##    Year State Class                 Strat_Category Data_Value
##   <dbl> <fct> <fct>                 <fct>               <dbl>
## 1  2019 Ohio  Fruits and Vegetables Race/Ethnicity       45.1

Zobaczmy jak sytuacja wygląda dla 10 rekordów:

dane_state %>% 
  mutate(CI_Width = High_Confidence_Limit -Low_Confidence_Limit) %>% 
  select(Year, State, Class, Strat_Category, CI_Width) %>% 
  arrange(desc(CI_Width)) %>% 
  head(10)

## # A tibble: 10 × 5
##     Year State    Class                   Strat_Category CI_Width
##    <dbl> <fct>    <fct>                   <fct>             <dbl>
##  1  2019 Ohio     Fruits and Vegetables   Race/Ethnicity     59.2
##  2  2019 Ohio     Obesity / Weight Status Race/Ethnicity     58.3
##  3  2023 Ohio     Obesity / Weight Status Race/Ethnicity     57.6
##  4  2023 Ohio     Physical Activity       Race/Ethnicity     55.7
##  5  2022 Virginia Physical Activity       Race/Ethnicity     55  
##  6  2015 Maryland Physical Activity       Race/Ethnicity     55  
##  7  2023 Ohio     Physical Activity       Race/Ethnicity     54.4
##  8  2023 Ohio     Obesity / Weight Status Race/Ethnicity     52.4
##  9  2023 Georgia  Physical Activity       Race/Ethnicity     52.2
## 10  2023 Ohio     Physical Activity       Race/Ethnicity     52.1

Jednak więcej rekordów niż przypuszczałem może mieć dość szeroki przedział ufności. Poszerzmy pole widzenia przyjmując, że powyżej 20 pkt procentowych dane mogą być już uważane za mało wiarygodne.

dane_state <- dane_state %>% 
  mutate(CI_Width = High_Confidence_Limit -Low_Confidence_Limit)

dane_state %>% 
  filter(CI_Width > 20) %>% 
  select(Year, State, Class, Strat_Category, CI_Width)

## # A tibble: 7,530 × 5
##     Year State    Class                   Strat_Category CI_Width
##    <dbl> <fct>    <fct>                   <fct>             <dbl>
##  1  2011 Alabama  Obesity / Weight Status Race/Ethnicity     32.5
##  2  2011 Alabama  Obesity / Weight Status Race/Ethnicity     33.3
##  3  2011 Alabama  Obesity / Weight Status Race/Ethnicity     25.6
##  4  2011 Alaska   Obesity / Weight Status Race/Ethnicity     24.7
##  5  2011 Alaska   Obesity / Weight Status Race/Ethnicity     28.2
##  6  2011 Arizona  Obesity / Weight Status Race/Ethnicity     26.9
##  7  2011 Arizona  Obesity / Weight Status Race/Ethnicity     23.3
##  8  2011 Arizona  Obesity / Weight Status Race/Ethnicity     23.7
##  9  2011 Arkansas Obesity / Weight Status Race/Ethnicity     29.6
## 10  2011 Arkansas Obesity / Weight Status Race/Ethnicity     23  
## # ℹ 7,520 more rows

Około 7.5 tysiąca rekordów. Ciekawe czy ponownie wpływ na to miałby grupy stratyfikacji czy coś innego. Rozbijmy to na przypadki by mieć szerszy pogląd sytuacji:

Kategoria

kat_wartosc <- dane_state %>% 
  filter(CI_Width > 20) %>% 
  count(Strat_Category, sort = TRUE )

kat_wartosc

## # A tibble: 4 × 2
##   Strat_Category     n
##   <fct>          <int>
## 1 Race/Ethnicity  7062
## 2 Income           301
## 3 Education        123
## 4 Age (years)       44

Ponad 7 tys. rekordów dotyczy pochodzenia. Pozostałe kategorie stanowią zdecydowaną mniejszość

Pytanie

pyt_wartosc <- dane_state %>% 
  filter(CI_Width > 20) %>% 
  count(Question, sort = TRUE )

pyt_wartosc

## # A tibble: 9 × 2
##   Question                                                                     n
##   <chr>                                                                    <int>
## 1 Percent of adults aged 18 years and older who have an overweight classi…  1526
## 2 Percent of adults aged 18 years and older who have obesity                1378
## 3 Percent of adults who engage in no leisure-time physical activity         1285
## 4 Percent of adults who achieve at least 150 minutes a week of moderate-i…   771
## 5 Percent of adults who engage in muscle-strengthening activities on 2 or…   731
## 6 Percent of adults who achieve more than 300 minutes a week of moderate-…   671
## 7 Percent of adults who achieve at least 150 minutes a week of moderate-i…   564
## 8 Percent of adults who report consuming fruit less than one time daily      355
## 9 Percent of adults who report consuming vegetables less than one time da…   249

Widać najbardziej problemowe pytania. Czy są one związane z kategorią?

Połączmy ostatnie dwie zmienne

zestawienie_kat_qest <- dane_state %>%
  filter(CI_Width > 20) %>%
  count(Question, Strat_Category, sort = TRUE)

zestawienie_kat_qest

## # A tibble: 36 × 3
##    Question                                                 Strat_Category     n
##    <chr>                                                    <fct>          <int>
##  1 Percent of adults aged 18 years and older who have an o… Race/Ethnicity  1460
##  2 Percent of adults aged 18 years and older who have obes… Race/Ethnicity  1292
##  3 Percent of adults who engage in no leisure-time physica… Race/Ethnicity  1205
##  4 Percent of adults who achieve at least 150 minutes a we… Race/Ethnicity   716
##  5 Percent of adults who engage in muscle-strengthening ac… Race/Ethnicity   686
##  6 Percent of adults who achieve more than 300 minutes a w… Race/Ethnicity   623
##  7 Percent of adults who achieve at least 150 minutes a we… Race/Ethnicity   533
##  8 Percent of adults who report consuming fruit less than … Race/Ethnicity   323
##  9 Percent of adults who report consuming vegetables less … Race/Ethnicity   224
## 10 Percent of adults aged 18 years and older who have obes… Income            59
## # ℹ 26 more rows

Czyli to głównie kategoria Race/Ethnicity ma mało wiarygodne dane.
Szybki podgląd na wykres:

zestawienie<- dane_state %>%
  filter(CI_Width > 20) %>%
  count(Question, Strat_Category) %>% 
  left_join(legenda_pytan, by = "Question") %>% 
  complete(ID, Strat_Category, fill = list(n = 0))

  zestawienie %>%
    group_by(Strat_Category) %>%
    summarise(
      n = sum(n), 
      .groups = "drop") %>%
  
    ggplot(aes(x=reorder(Strat_Category, -n), y = n)) +
    geom_col(fill = "steelblue") +
    labs(
      title = "Liczba przypadków z CI > 20 według kategorii",
      x = "Kategoria stratyfikacji",
      y = "Liczba przypadków"
    ) +
    theme_minimal() +
    theme(
      text = element_text(size = 12),
      plot.title = element_text(size = 15, face = "bold"),
      axis.title = element_text(size = 14),
    )

Patrząc na problemy związane z Race/Ethnicity zasadne wydaje się pominięcie tej zmiennej w dalszej analizie.

Dodatkowo przeanalizujmy rozkład danych w ujęciu czasowym:

rok_wartosc <- dane_state %>% 
  filter(CI_Width > 20) %>% 
  count(Year) %>% 
  
  ggplot(aes(x = as.factor(Year), y = n)) +
  geom_col(fill = "steelblue") +
  labs(
    title = "Liczba niepewnych rekordów na przestrzeni lat",
    subtitle = "CI Width > 20",
    x = "Rok",
    y = "Liczba wystąpień"
  ) +
  theme_bw()+
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

rok_wartosc

Różna liczba rekordów o dużym przedziale ufności, może wpływać na wiarygodność porównań między latami czy po prostu nawet podczas wyliczania średniej. Jednak nie zapominajmy, że nie wszystkie pytania były zadawane co roku i to jest główną przyczyną dlaczego w latach nieparzystych jest więcej rekordów.

Przyjrzyjmy się jeszcze naszym wynikom (Data_Value)

summary(dane_state$Data_Value)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.90   24.80   31.70   31.76   37.40   85.30   10226

Średnia ma prawie tę samą wartość co mediana - rozkład jest prawie symetryczny (wartość maksymalna jest całkiem mocno oddalona od 3 kwartyla, więc to ona mogła pociągnąć średnią do góry)
Rozstęp międzykwartylowy IQR = 37.4 - 31.7 = 5.7 - co świadczy o niskiej zmienności w obrębie środkowych 50% obserwacji.
Interesuje mnie ta wartość maksymalna - 85.30. Zobaczmy na jakie pytanie ludzie tak licznie odpowiedzieli twierdząco

dane_state %>%
  filter(Data_Value > 70) %>%
  select(Question, State, Data_Value, Strat_Group) %>%
  arrange(desc(Data_Value)) %>%
  head(10) %>%
  knitr::kable(caption = "Obserwacje z odsetkiem > 70%" )

Obserwacje z odsetkiem > 70%
Question	State	Data_Value	Strat_Group
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	Texas	85.3	American Indian/Alaska Native
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	District of Columbia	80.0	$75,000 or greater
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	District of Columbia	79.3	Asian
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	Oregon	77.6	Asian
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	District of Columbia	76.8	Non-Hispanic White
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	Vermont	76.6	College graduate
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	Maine	76.3	College graduate
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	Alaska	76.2	College graduate
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	Idaho	75.5	College graduate
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	Montana	75.3	Other

To pytanie dotyczy osetka dorosłych, którzy spełniają minimalne rekomendacje WHO dotyczące aktywności fizycznej (≥150 min umiarkowanej lub ≥75 min intensywnej aktywności tygodniowo). Nie jest to wynik dziwny, ponieważ nawet zwykłe codzienne spacery mogą się przyczynić do wyrobienia zalecanej normy.

Oczywiście rozkład zmiennej Data_Value wynika przede wszystkim z różnorodności pytań i tak naprawdę ciężko mówić o wartościach odstających, skoro dana wartość może dotyczyć różnego pytania. Lepiej więc patrzeć na dane pytanie z osobna lub badać stopień zróżnicowania (np między ilością aktywności fizycznej)

Mając już wiedzę o strukturze zbioru, przejdźmy do odpowiedzi na kluczowe pytania badawcze…

W którym rejonie odsetek osób z otyłością jest największy?

Czy faktycznie zachód cechuje się najniższymi wskaźnikami otyłości, podczas gdy wschód posiada najgorsze wskaźniki w kraju?

dane1 <- dane_state %>%
  filter(Question == "Percent of adults aged 18 years and older who have obesity") %>% 
  filter(Strat_Category == "Total") %>%
  group_by(Region) %>%
  summarise(
    suma = sum(Data_Value, na.rm = TRUE),
    liczba = sum(!is.na(Data_Value)),
    srednia = suma / liczba
  ) 

knitr::kable(
  dane1 %>%  select(Region, srednia), 
  caption = "Średni odsetek otyłości wg regionu")

Średni odsetek otyłości wg regionu
Region	srednia
North Central	32.78155
Northeast	27.98871
South	33.35404
West	28.04560

dane1 %>% 
  ggplot(aes(x=reorder(Region, -srednia), y = srednia)) +
  geom_col(fill = "steelblue") +
  labs(
    title = "Średni odsetek otyłości wg regionu",
    x = "Region",
    y = "Średni oodsetek (%)"
  ) +
  theme_bw()+
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

Widzimy, że nie ma dużych różnic pomiędzy regionami.
Najmniejszym współczynnikiem cieszy się Northeast = ~28%, za to największym South = 33.35%

Dla potwierdzenia zastosujmy test statystyczny:
Ho: Regiony nie różnią się poziomem otyłości
H1: Przynajmniej jeden region ma inny poziom otyłości

kw <- dane_state %>%
  filter(Question == "Percent of adults aged 18 years and older who have obesity") %>% 
  filter(Strat_Category == "Total") %>% 
  group_by(State, Region) %>%
  summarise(
    obesity = mean(Data_Value), 
    .groups = "drop")

kruskal.test(obesity ~ Region, data = kw)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  obesity by Region
## Kruskal-Wallis chi-squared = 26.638, df = 3, p-value = 7.01e-06

Ponieważ p-value < 0.05 to odrzucamy hipotezę zerową.
Potwierdziliśmy, że co najmniej jeden region ma inny średni poziom otyłości.

Na nasze wyniki mogą mieć wpływ czynniki ekonomiczne i edukacyjne. Historycznie południe to region o najwyższym wskaźniku ubóstwa i najniższych dochodach. Przy czym odwrotnie na północy są najwyższe zarobki i lepszy poziom edukacji.

W takim razie zobaczmy jak konkretne czynniki wpływają na otyłość:

Czy niższy poziom dochodu wiąże się z wyższym odsetkiem osób z otyłością

Z jednej strony wyższe dochody mogą sprzyjać lepszym nawykom żywieniowym i większej dbałości o jakość diety. Z drugiej strony, praca biurowa - często skorelowana z wyższymi zarobkami – wiąże się z osiadłym trybem życia, co stanowi istotny czynnik ryzyka dla utrzymania prawidłowej masy ciała. Sprawdźmy, co na ten temat mówią zgromadzone dane.

Na początku przygotujmy poziomy dochodu

dochody<- c(
  "Data not reported",
  "Less than $15,000",
  "$15,000 - $24,999",
  "$25,000 - $34,999",
  "$35,000 - $49,999",
  "$50,000 - $74,999",
  "$75,000 or greater"
)

dane2 <- dane_state %>%
  filter(Strat_Category == "Income",
    Question == "Percent of adults aged 18 years and older who have obesity"
  ) %>%
  mutate(Strat_Group = factor(Strat_Group, levels = dochody )) %>% 
  group_by(Region, Strat_Group) %>% 
  summarise(
    srednia = mean(Data_Value, na.rm = TRUE),
    n = sum(!is.na(Data_Value)),
    .groups = "drop"
  )

dane2 %>% 
  ggplot(aes(x = srednia, y = Strat_Group)) +
  geom_segment(aes(x = 0, xend = srednia, yend = Strat_Group),
               linewidth = 0.6, alpha = 0.5) +
  geom_point(size = 3) +
  facet_wrap(~Region) +
  labs(
    title = "Poziom otyłości a dochód",
    x = "Średni % otyłości",
    y = "Przedział dochodowy"
  ) +
  theme_minimal()+  
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

Dane potwierdzają istnienie gradientu zdrowotnego: wraz ze wzrostem statusu materialnego, wskaźnik otyłości ulega obniżeniu. Zjawisko to można przypisać większej dostępności zasobów sprzyjających zdrowemu stylowi życia, takich jak wysokiej jakości produkty spożywcze, personalizowane usługi dietetyczne (catering dietetyczny) czy profesjonalne wsparcie w zakresie aktywności fizycznej.

Jaki jest związek poziomu wykształcenia z odsetkiem osób z otyłością?

Wpływ wykształcenia na wskaźnik otyłości można rozpatrywać dwutorowo. Z jednej strony, proces edukacji często wiąże się z wieloletnim, siedzącym trybem życia i deficytem aktywności fizycznej. Z drugiej strony, wyższe wykształcenie może wpływać na większą świadomością prozdrowotną, i tym samym zdrowszy styl życia.

wyksztalcenie <- c(
  "Less than high school",
  "High school graduate",
  "Some college or technical school",
  "College graduate"
)

dane3 <- dane_state %>%
  filter(Strat_Category == "Education") %>% 
  filter(Question == "Percent of adults aged 18 years and older who have obesity") %>% 
  mutate(Strat_Group = factor(Strat_Group, levels = wyksztalcenie)) %>% 
  group_by(Region, Strat_Group) %>% 
  summarise(
    srednia = mean(Data_Value, na.rm = TRUE),
    n = sum(!is.na(Data_Value)),
    .groups = "drop"
  )

dane3 %>% 
  ggplot(aes(x = srednia, y = Strat_Group)) + 
  geom_segment(aes(x = 0, xend = srednia, yend = Strat_Group),
               linewidth = 0.6, alpha = 0.5) +
  geom_point(size = 3) +
  facet_wrap(~Region) +
  labs(
    title = "Poziom otyłości a wykstałcenie",
    x = "Średni % otyłości",
    y = "Rodzaj wykształcenia"
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

Widzimy, że ludzie, którzy nie ukończyli college’u, cechują się wyższym poziomem otyłości. Dla pozostałych 3 kategorii odsetek oscyluje w granicach 30-37%. Wcześniej wykazaliśmy, że wyższy dochód jest czynnikiem wpływającym na otyłość, a zazwyczaj ludzie z wyższym wykształceniem zarabiają więcej.

Wspomnieliśmy wcześniej o siedzącym trybie życia. Jest powszechnie wiadomo, że im więcej się ruszasz tym więcej kalorii spalasz. Czyli względnie ludzie, którzy ruszają się więcej mają mniejsze szanse na zostanie otyłym. Choć oprócz samej ilości znacząca jest także intensywność:

W ciągu 15 minut:

w wolnym biegu spalimy ok 130-150 kcal
Podczas skakania na skakance ok 190-220 kcal

W takim razie przyjrzyjmy się tej kwestii bliżej: Na początku zobaczmy jak ogólnie Amerykanie podchodzą do aktywności fizycznej, a następnie sprawdźmy jak sama aktywność wiąże się z otyłością.

Jaki rodzaj aktywności fizycznej jest najczęściej deklarowany przez dorosłych?

Czy większość społeczeństwa spełnia minimalne normy aktywności fizycznej WHO?

Na początku zobaczmy jaki mamy podział aktywności fizycznej

aktywnosc <- dane_state %>% 
  filter(Class == "Physical Activity") %>% 
  distinct(Question)

aktywnosc

## # A tibble: 5 × 1
##   Question                                                                      
##   <chr>                                                                         
## 1 Percent of adults who achieve at least 150 minutes a week of moderate-intensi…
## 2 Percent of adults who achieve at least 150 minutes a week of moderate-intensi…
## 3 Percent of adults who achieve more than 300 minutes a week of moderate-intens…
## 4 Percent of adults who engage in muscle-strengthening activities on 2 or more …
## 5 Percent of adults who engage in no leisure-time physical activity

Widać że tutaj również warto przypisać poziomy

poz_aktywnosc <- c(
  "Percent of adults who engage in no leisure-time physical activity",
  "Percent of adults who engage in muscle-strengthening activities on 2 or more days a week",
  "Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)",
  "Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic physical activity (or an equivalent combination) and engage in muscle-strengthening activities on 2 or more days a week",
  "Percent of adults who achieve more than 300 minutes a week of moderate-intensity aerobic physical activity or 150 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)"
)

Z uwagi na czytelność wykresu, przypisze pytaniom (a właściwie im poziomom) skróconą nazwę

legenda_aktywnosci <- data.frame(Question = poz_aktywnosc,Level = 1:5) %>%
  mutate(Level = factor(Level))

knitr::kable(legenda_aktywnosci, caption = "Poziomy aktywności")

Poziomy aktywności
Question	Level
Percent of adults who engage in no leisure-time physical activity	1
Percent of adults who engage in muscle-strengthening activities on 2 or more days a week	2
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	3
Percent of adults who achieve at least 150 minutes a week of moderate-intensity aerobic physical activity or 75 minutes a week of vigorous-intensity aerobic physical activity (or an equivalent combination) and engage in muscle-strengthening activities on 2 or more days a week	4
Percent of adults who achieve more than 300 minutes a week of moderate-intensity aerobic physical activity or 150 minutes a week of vigorous-intensity aerobic activity (or an equivalent combination)	5

Myślę, że warto wytłumaczyć skąd taka hierarchia:

Brak aktywności: najwyższe ryzyko przyrostu masy ciała i chorób sercowo-naczyniowych.
Tylko ćwiczenia siłowe: wzmacniasz mięśnie i kości, poprawiasz sylwetkę i metabolizm, ale mniej rozwijasz wydolność tlenową.
Minimum cardio (zalecenia WHO): spełniasz podstawowe normy zdrowotne: 150–300 min/tydz umiarkowanej aktywności (szybki marsz) lub 75–150 min/tydz intensywnej (bieg)
Cardio + siłowe: najbardziej “kompletna” opcja – łączy korzyści dla wydolności (serce, układ krążenia) i dla mięśni/kości.
Powyżej 300 min/tydz.: bardzo wysoka aktywność (budowanie kondycji, znaczne spalanie kalorii)

Z uwagi na dostępność danych do analizy weźmiemy dane z lat: 2011, 2013, 2015, 2017, 2019, 2023

Najpierw zobaczmy jak to wygląda zbiorczo

badane_lata <- c(2011, 2013, 2015, 2017, 2019, 2023)

dane4 <- dane_state %>%
  filter(Year %in% badane_lata) %>% 
  filter(Class == "Physical Activity") %>% 
  filter(Strat_Category == "Total") %>% 
  group_by(Question, Region) %>% 
    summarise(
    suma = sum(Data_Value, na.rm = TRUE),
    liczba = sum(!is.na(Data_Value)),
    srednia = suma / liczba,
    .groups = "drop"
  )

dane4 <- dane4 %>%
  left_join(legenda_aktywnosci, by = "Question")

dane4 %>% 
  ggplot() +
  geom_col(aes(x = Region, y = srednia, fill = Level),
           color = "white",
           position = "dodge2") +
  scale_fill_brewer(palette = "Blues")  +
  theme_bw() +
  labs(title = "Poziom aktywności fizycznej w poszczególncyh regionach",
       x = "Region",
       y = "Średnia wartość (%)",
       fill = "Poziom (1-5)") +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

Spostrzeżenia: W każdym rejonie balans poziomu aktywności jest podobny:

Średnie dla poziomu 3 wypadają najlepiej - duży odsetek ludzi spełnia minimum zalecane przez WHO.
Najmniejszy odsetek dotyczy najbardziej kompletnej opcji - cardio + siłowe (jedynie wyjątek w rejonie West, gdzie najmniejszy odsetek dotyczy braku poświęcania wolnego czasu na aktywność fizyczną)
W każdym rejonie ponad 20% społeczeństwa nie poświęca wolnego czasu na aktywność fizyczną - jest to dość sporo.
Aż 30-40% wykazuje intensywną aktywność fizyczną (level 5)

Czy wyższy poziom intensywności aktywności fizycznej wiąże się z niższym odsetkiem otyłości?

Sprawdźmy teraz czy powszechna wiedza: więcej ruchu = lepsza sylwetka, ma potwierdzenie w danych

badane_lata <- c(2011, 2013, 2015, 2017, 2019, 2023)

dane_otylosc <- dane_state %>%
  filter(Year %in% badane_lata) %>% 
  filter(Question == "Percent of adults aged 18 years and older who have obesity") %>% 
  filter(Strat_Category == "Total") %>% 
  select(Year, State, otylosc_proc = Data_Value)

dane_sport <- dane_state %>%
  filter(Year %in% badane_lata) %>% 
  filter(Class == "Physical Activity") %>% 
  filter(Strat_Category == "Total") %>% 
  select(Year, State, Question, aktywnosc_proc = Data_Value) %>% 
  left_join(legenda_aktywnosci, by = "Question")

dane_sport_szerokie <- dane_sport %>%
  select(Year, State, Level, aktywnosc_proc) %>% 
  pivot_wider(
    names_from = Level,
    names_prefix = "L",
    values_from = aktywnosc_proc
  )

dane5 <- dane_sport_szerokie %>% 
  left_join(dane_otylosc, by = c("Year", "State"))

Korelacja Spearmana

korelacja_sp <- dane5 %>% 
  select("L1", "L2", "L3", "L4", "L5") %>% 
  sapply(cor, y = dane5$otylosc_proc, method = "spearman", use = "complete.obs") %>% 
  sort() 

knitr::kable(korelacja_sp, col.names = "Korelacja", caption = "Związek aktywności z otyłością")

Związek aktywności z otyłością
	Korelacja
L3	-0.3385203
L4	-0.2335127
L5	-0.2231121
L2	-0.0954523
L1	0.5586948

Zaobserwowaliśmy dodatnią zależność dla osób niepoświęcających czasu na aktywność fizyczną oraz coraz bardziej ujemną zależność wraz ze zwiększeniem stopnia aktywności fizycznej. Czyli wyższy poziom aktywności jest powiązany z niższym odsetkiem otyłości.

Wcześniej badaliśmy jak wykształcenie wpływa na poziom otyłości. (dla przypomnienia w skrócie: im wyższy poziom wykształcenia, tym mniejsza otyłość). Skoro przed chwilą zobaczyliśmy wpływ ruchu na poziom naszej tkanki tłuszczowej to zobaczmy czy niższy poziom otyłości wynika z większej ilości ruchu.

Czy poziom wykształcenia jest związany z poziomem aktywności fizycznej?

dane7 <- dane_state %>% 
  filter(Year %in% badane_lata) %>% 
  filter(Class == "Physical Activity") %>% 
  filter(Strat_Category == "Education") %>% 
  select(Year, State, Question, Data_Value, Strat_Group) %>% 
  left_join(legenda_aktywnosci, by = "Question")

dane7.1 <- dane7 %>%
  mutate(Strat_Group = factor(Strat_Group, levels = wyksztalcenie)) %>%
  group_by(Strat_Group, Level) %>% 
    summarise(
    suma = sum(Data_Value, na.rm = TRUE),
    liczba = sum(!is.na(Data_Value)),
    srednia = suma / liczba,
    .groups = "drop"
    )

dane7.1 %>% 
  ggplot() +
  geom_col(aes(x= Strat_Group, y = srednia, fill = Level),
           color = "white",
           position = "dodge2") +
  scale_fill_brewer(palette = "Blues")  +
  theme_bw() +
  labs(
    title = "Poziom wykształcenia, a ilość aktywności fizycznej",
    x = "Wykształcenie",
    y = "Odsetek badanych [%]",
    fill = "Poziom aktywności"
  )+
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

Wraz ze wzrostem poziomu wykształcenia obserwujemy wyraźny wzrost ogólnej aktywności fizycznej. Odsetek osób całkowicie unikających sportu w grupie z wykształceniem wyższym jest ponad dwukrotnie niższy niż wśród osób bez wykształcenia. Jednocześnie wraz ze wzrostem poziomu wykształcenia wzrasta odsetek osób deklarujący intensywny poziom aktywności (Poziom 4 oraz 5)

W naszym zbiorze danych mamy również dane dotyczące diety:

procent jedzący mniej niej 1 owoc dziennie
procent jedzący mniej niż 1 warzywo dziennie

Szybkie rozeznanie:

owoce <- dane_state %>% 
  filter(Class == "Fruits and Vegetables") %>% 
  filter(Question == "Percent of adults who report consuming fruit less than one time daily") %>% 
  filter(Strat_Category == "Total") %>%
  select(Data_Value)
  
summary(owoce)

##    Data_Value   
##  Min.   :29.70  
##  1st Qu.:36.05  
##  Median :38.30  
##  Mean   :39.01  
##  3rd Qu.:42.55  
##  Max.   :48.70  
##  NA's   :2

warzywa <- dane_state %>% 
  filter(Class == "Fruits and Vegetables") %>% 
  filter(Question == "Percent of adults who report consuming vegetables less than one time daily") %>% 
  filter(Strat_Category == "Total") %>%
  select(Data_Value)
  
summary(warzywa)

##    Data_Value   
##  Min.   :12.40  
##  1st Qu.:17.25  
##  Median :19.40  
##  Mean   :19.29  
##  3rd Qu.:21.30  
##  Max.   :25.70  
##  NA's   :2

pytania_dieta <- c("Percent of adults who report consuming fruit less than one time daily","Percent of adults who report consuming vegetables less than one time daily")

warzywaIowoce <- dane_state %>%
  filter(Strat_Category == "Total") %>% 
  filter(Question %in% pytania_dieta)

warzywaIowoce %>% 
  ggplot(aes(x = Question, y = Data_Value, fill = Question)) +
  geom_boxplot(width = 0.6, alpha = 0.7) +
  scale_x_discrete(labels = c(
    "Percent of adults who report consuming fruit less than one time daily" = "Owoce",
    "Percent of adults who report consuming vegetables less than one time daily" = "Warzywa"
  )) +
  scale_fill_brewer(palette = "Set2") +
  labs(
    title = "Porównanie niskiego spożycia owoców i warzyw",
    x = "",
    y = "Odsetek (%)"
  ) +
  theme_bw() +
  theme(legend.position = "none")+
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

Widać wyraźną różnicę między badanymi kategoriami: odsetek osób spożywających owoce rzadziej niż raz dziennie jest istotnie wyższy niż w przypadku warzyw. Mediana dla owoców wynosi około 38–39%, podczas gdy dla warzyw około 18–19%, co oznacza niemal dwukrotnie wyższy poziom niskiego spożycia owoców. Dodatkowo rozkład dla owoców charakteryzuje się większą zmiennością, co sugeruje silniejsze zróżnicowanie regionalne w porównaniu do warzyw.

Przejdźmy do naszego pytania

Który czynnik jest silniej związany z otyłością: brak aktywności fizycznej czy niskie spożycie owoców?

Na początku zobaczmy jaki wpływ na otyłość mają same owoce:

(Pytania o owoce i warzywa były zadawane w latach: 2017, 2019, 2021)

badane_lata_owoce <- c(2017, 2019, 2021)

dane_otylosc_6 <- dane_state %>%
  filter(Year %in% badane_lata_owoce) %>% 
  filter(Question == "Percent of adults aged 18 years and older who have obesity") %>% 
  filter(Strat_Category == "Total") %>% 
  select(Year, State, otylosc_proc = Data_Value)

dane_owoce <- dane_state %>%
  filter(Year %in% badane_lata_owoce) %>% 
  filter(Class == "Fruits and Vegetables") %>% 
  filter(Question == "Percent of adults who report consuming fruit less than one time daily") %>% 
  filter(Strat_Category == "Total") %>% 
  select(Year, State, Question, owoc_procent = Data_Value)

dane6 <- dane_otylosc_6 %>% 
  left_join(dane_owoce, by = c("Year", "State"))

korelacja_ruch <- dane_state %>% 
  filter(Year %in% badane_lata_owoce) %>% 
  filter(Question == "Percent of adults who engage in no leisure-time physical activity") %>% 
  filter(Strat_Category == "Total") %>% 
  select(Year, State, Question, bezruch_procent = Data_Value)

dane6.1 <- dane_otylosc_6 %>% 
  left_join(korelacja_ruch, by = c("Year", "State"))

ggplot(dane6, aes(x = owoc_procent, y = otylosc_proc)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  facet_wrap(~ Year) +
  labs(
    title = "Spożycie owoców a otyłość",
    x = "Odsetek jedzących owoce rzadziej niż 1× dziennie",
    y = "Otyłość (%)"
  ) +
  theme_bw()

cor(dane6$owoc_procent, dane6$otylosc_proc, use = "complete.obs", method = "spearman")

## [1] 0.7723239

cor(dane6.1$bezruch_procent, dane6$otylosc_proc, use = "complete.obs", method = "spearman")

## [1] 0.5749459

Potwierdźmy nasze wyniki testem:
Ho: p = 0 (brak zależności)
H1: p ≠ 0 (istnieje zależność)

cor.test(dane6.1$otylosc_proc, dane6.1$bezruch_procent,
         method = "spearman", use = "complete.obs")

## 
##  Spearman's rank correlation rho
## 
## data:  dane6.1$otylosc_proc and dane6.1$bezruch_procent
## S = 243896, p-value = 1.161e-14
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.5749459

Odrzucamy hipotezę zerową Ho.
W danych występuje istotna statystycznie dodatnia korelacja.

Zarówno wskaźnik braku aktywności, jak i rzadszego spożycia owoców są dodatnio skorelowane z otyłością. Jednak korelacja otyłości z rzadszym spożyciem owoców jest silniejsza niż z brakiem aktywności w czasie wolnym. Z perspektywy praktycznej łatwiej jest zwiększyć podaż energii (np. o 300 kcal) niż ją spalić.

Skoro już jesteśmy przy temacie odżywiania i wcześniej sprawdzaliśmy, jaki wpływ ma wykształcenie na różne czynniki to naturalnym następnym krokiem będzie sprawdzenie jak poziom wykształcenia wpływa na nasze nawyki żywieniowe.

Czy poziom dochodu jest związany z częstotliwością spożycia owoców i warzyw?

dane9 <- dane_state %>% 
  filter(Question == "Percent of adults who report consuming fruit less than one time daily" | Question == "Percent of adults who report consuming vegetables less than one time daily") %>% 
  filter(Strat_Category == "Income") %>% 
  mutate(Strat_Group = factor(Strat_Group, levels = dochody )) %>% 
  group_by(Question, Strat_Group) %>% 
  summarise(
    srednia = mean(Data_Value, na.rm = TRUE),
    .groups = "drop"
  )

dane9 %>% 
  ggplot(aes(x = Strat_Group, y = srednia))+
  geom_col(fill = "steelblue") +
  coord_flip() +
  facet_wrap(~Question, ncol = 1) +
  labs(
    title = "Poziom dochodu a nawyki żywieniowe",
    x= "Poziom dochodu",
    y= "Odsetek osób jedzących < 1 raz dziennie"
  ) + 
  theme_bw()+
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

Dzięki wizualizacji widzimy, że wyższy poziom dochodu sprzyja zdrowszym decyzjom żywieniowym. Jednak większa różnica jest w konsumpcji warzyw - około 16.65%, podczas gdy różnica w konsumpcji owoców wynosi około 8,6%

Czy w regionach otyłość rośnie, maleje czy jest stabilna?

dane8 <- dane_state %>% 
  filter(Question == "Percent of adults aged 18 years and older who have obesity") %>% 
  filter(Strat_Category == "Total") %>%
  group_by(Year, Region) %>% 
  summarise(
    srednia = mean(Data_Value, na.rm = TRUE),
    .groups = "drop"
  )

dane8 %>% 
  ggplot()+
  geom_line(aes(x = as.factor(Year), y = srednia, color = Region, group = Region),
            linewidth = 1.2) +
  scale_color_brewer(palette = "Set1")+
  labs(
    title = "Poziom otyłości z biegiem lat",
    x= "Lata",
    y= "Średni poziom otyłości [%]",
    color= "Region"
  )+
  theme_bw() +
  theme(
    text = element_text(size = 12),
    plot.title = element_text(size = 15, face = "bold"),
    axis.title = element_text(size = 14),
  )

W latach 2011–2024 obserwujemy wyraźny wzrost odsetka dorosłych osób z otyłością we wszystkich regionach USA. Najwyższe poziomy otyłości przez cały okres charakteryzują regiony South oraz North Central, które wyraźnie odstają od pozostałych części kraju. Warto jednak zauważyć, że region South jako jedyny wykazuje w ostatnich dwóch latach niewielką tendencję spadkową, podczas gdy w pozostałych regionach poziom otyłości utrzymuje tendencje wzrostową.