Body signal of smoking

Wstęp

    Jako projekt wybrałem analizę danych, który jest zbiorem podstawowych informacji dotyczących biologicznych sygnałów zdrowotnych. Celem jest określenie czy określona osoba jest palaczem czy nie za pomocą biosygnałów.

Elementy projektu:

• opis danych,
• przygotowanie danych,
• analiza zbioru danych
• PCA,
• Drzewa decyzyjne,
• Neural networks,
  • original data,
  • PCA data,
• Wnioski i podsumowanie.

     Poniżej zaprezentowane są dane na których będę pracował. Zbiór zawiera 26 zmiennych na temat poszczególnych osób, palących i niepalących.

df <- read.csv2("smoking.csv", sep = ",", header = TRUE)
head(df, 2)

##   ID gender age height.cm. weight.kg. waist.cm. eyesight.left. eyesight.right.
## 1  0      F  40        155         60      81.3            1.2             1.0
## 2  1      F  40        160         60      81.0            0.8             0.6
##   hearing.left. hearing.right. systolic relaxation fasting.blood.sugar
## 1           1.0            1.0    114.0       73.0                94.0
## 2           1.0            1.0    119.0       70.0               130.0
##   Cholesterol triglyceride  HDL   LDL hemoglobin Urine.protein serum.creatinine
## 1       215.0         82.0 73.0 126.0       12.9           1.0              0.7
## 2       192.0        115.0 42.0 127.0       12.7           1.0              0.6
##    AST  ALT  Gtp oral dental.caries tartar smoking
## 1 18.0 19.0 27.0    Y             0      Y       0
## 2 22.0 19.0 18.0    Y             0      Y       0

Możemy zobaczyć kolumny:

• ID : index
• gender
• age : zaokrąglone do 5 letnich odstępów (5, 10, 15…)
• height(cm)
• weight(kg)
• waist(cm) : długość obwodu talii
• eyesight(left): wzrok lewego oka
• eyesight(right): wzrok prawego oka
• hearing(left): słuch, lewe ucho
• hearing(right): słuch, prawe ucho
• systolic : ciśnienie skurczowe
• relaxation : ciśnienie krwi
• fasting blood sugar: cukier we krwi na czczo
• Cholesterol:
• triglyceride
• HDL : typ cholesterolu
• LDL : typ cholesterolu
• hemoglobin
• Urine protein
• serum creatinine: białko w moczu
• AST : typ transaminazy glutaminowo-szczawiooctowej
• ALT : typ transaminazy glutaminowo-szczawiooctowej
• Gtp : γ-GTP
• oral :
• dental caries
• tartar : stan kamienia nazębnego
• smoking

Przygotowanie danych

Ładowanie potrzebnych bibliotek:

library(tree)
require(caret)
library(nnet)
library(corrplot)
library(reshape2)
library(dplyr)
set.seed(21)

Usuwanie pustych wierszy:

df <- df[complete.cases(df), ]   

Zamiana kolumny gender, oral oraz tartar z warości tekstowych na liczbowe (0, 1):

unique_gender <- unique(as.character(df$gender))
df <- data.frame(gender = as.numeric(factor(df$gender, level = unique_gender)), df)

unique_oral <- unique(as.character(df$oral))
df <- data.frame(oral = as.numeric(factor(df$oral, level = unique_oral)), df)

unique_tartar <- unique(as.character(df$tartar))
df <- data.frame(tartar = as.numeric(factor(df$tartar, level = unique_tartar)), df)

Usunięcie zbędnych zmiennych:

df <- subset(df, select = -c(ID, oral.1, gender.1, tartar.1))

Zamiana typu kolumn z character na numeric:

df[] <- lapply(df, function(x) as.numeric(as.character(x)))
sapply(df, class)

Analiza zbioru danych

    W pierwszej kolejności została zbadana korelacja pomiędzy badanymi zmiennymi:

corrplot(cor(df[,c(-2)], method='spearman'),
  method = "number",
  type = "upper")

    Jak widać zmienne nie są pomiędzy sobą mocno skorelowane. Można dostrzec, że płeć jest najbardziej skorelowana z paleniem, ale także jest najbardziej powiązana z resztą zmiennych. Po za tym palenie najmocniej skorelowane jest z hemoglobiną (0,40) oraz wzrostem (0,40).

    Następnie zostały przedstawione zmienne dyskretne jako wykresy kołowe:

p1 <-df%>%
  count(tartar)%>%
  ggplot(aes(x = "", y = n ,fill = factor(tartar))) +
  geom_col() +
  coord_polar(theta = "y")+
  ggtitle("Rozkład zmiennej tartar") + 
  theme(plot.title = element_text(hjust = 0.5))

p2 <- df%>%
  count(oral)%>%
  ggplot(aes(x = "", y = n ,fill = factor(oral))) +
  geom_col() +
  coord_polar(theta = "y")+
  ggtitle("Rozkład zmiennej oral") + 
  theme(plot.title = element_text(hjust = 0.5))

p3 <- df%>%
  count(gender)%>%
  ggplot(aes(x = "", y = n ,fill = factor(gender))) +
  geom_col() +
  coord_polar(theta = "y")+
  ggtitle("Rozkład zmiennej gender") + 
  theme(plot.title = element_text(hjust = 0.5))

p4 <- df%>%
  count(smoking)%>%
  ggplot(aes(x = "", y = n ,fill = factor(smoking))) +
  geom_col() +
  coord_polar(theta = "y")+
  ggtitle("Rozkład zmiennej smoking") + 
  theme(plot.title = element_text(hjust = 0.5))

p5 <- df%>%
  count(dental.caries)%>%
  ggplot(aes(x = "", y = n ,fill = factor(dental.caries))) +
  geom_col() +
  coord_polar(theta = "y")+
  ggtitle("Rozkład zmiennej dental.caries") + 
  theme(plot.title = element_text(hjust = 0.5))

    Można z nich wywnioskować, że zmienna oral (ocena stanu ustnego) jest bezużyteczna w naszym zbiorze danych i możemy ją usunąć:

df <- subset(df, select = -c(oral))

    Następnie zostały przedstawione pozostałe zmienne jako rozkłady dla osób palących i niepalących:

p5 <- df %>%
  ggplot( aes(x = factor(smoking), y = age , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej wiek dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p6 <- df %>%
  ggplot( aes(x = factor(smoking), y = height.cm. , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej height.cm. dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p7 <- df %>%
  ggplot( aes(x = factor(smoking), y = weight.kg. , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej weight.kg. dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p8 <- df %>%
  ggplot( aes(x = factor(smoking), y = waist.cm. , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej waist.cm. dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p9 <- df %>%
  ggplot( aes(x = factor(smoking), y = systolic , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej systolic dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p10 <- df %>%
  ggplot( aes(x = factor(smoking), y = relaxation , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej relaxation dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p11 <- df %>%
  ggplot( aes(x = factor(smoking), y = fasting.blood.sugar , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej fasting.blood.sugar dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p12 <- df %>%
  ggplot( aes(x = factor(smoking), y = Cholesterol , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej Cholesterol dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p13 <- df %>%
  ggplot( aes(x = factor(smoking), y = triglyceride , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej triglyceride dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p14 <- df %>%
  ggplot( aes(x = factor(smoking), y = HDL , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej HDL dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p15 <- df %>%
  ggplot( aes(x = factor(smoking), y = LDL , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej LDL dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p16 <- df %>%
  ggplot( aes(x = factor(smoking), y = hemoglobin , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej hemoglobin dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p17 <- df %>%
  ggplot( aes(x = factor(smoking), y =AST  , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej AST dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p18 <- df %>%
  ggplot( aes(x = factor(smoking), y = ALT , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej ALT dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

p19 <- df %>%
  ggplot( aes(x = factor(smoking), y = Gtp , fill=factor(smoking))) +
    geom_violin() + 
    ggtitle("Rozkład zmiennej Gtp dla osób niepalący i palących ") +
    theme(plot.title = element_text(hjust = 0.5)) +
    ylab("")+
    xlab("") + 
    theme(legend.position="none")

    Z wykresów można odczytać, że w przypadku zmiennych: fasting.blood.sugar, HDL, LDL, AST, ALT, Gtp występują tak zwane outliery, które przyjmują wartości nierealne jeżeli chodzi o tego typu wyniki badań. Mogą one zaburzać wyniki, dlatego W celu normalizacji danych dla powyższych zmiennych zostaną usunięte dane, które są powyżej wybranego percentyla, a więc tylko dla bardzo małej liczby obserwacji.

attach(df)
q_fasting.blood.sugar <- quantile(fasting.blood.sugar, probs=c( .99), na.rm = FALSE)
q_HDL <- quantile(fasting.blood.sugar, probs=c( .99), na.rm = FALSE)
q_LDL <- quantile(fasting.blood.sugar, probs=c( .999), na.rm = FALSE)
q_AST <- quantile(fasting.blood.sugar, probs=c( .85), na.rm = FALSE)
q_ALT <- quantile(fasting.blood.sugar, probs=c( .90), na.rm = FALSE)
q_Gtp <- quantile(fasting.blood.sugar, probs=c( .99), na.rm = FALSE)

df <- subset(df, df$fasting.blood.sugar  < (q_fasting.blood.sugar[1]))
df <- subset(df, df$HDL < (q_HDL[1]))
df <- subset(df, df$LDL < (q_LDL[1]))
df <- subset(df, df$AST < (q_AST[1]))
df <- subset(df, df$ALT < (q_ALT[1]))
df <- subset(df, df$Gtp < (q_Gtp[1]))

    Po redukcji można zobaczyć jak teraz wyglądają rozkłądy poszczegółnych zmiennych:

    Jak widać po redukcji danych mocno odstających rozkłady zmiennych zmieniły się na bardziej skoncentrowane.

    Po analizie zbioru danych możesz przejść do podzielenia danych na train i test:

idx <- floor(0.75 * nrow(df))
train_ind <- sample(seq_len(nrow(df)), size = idx)

trainset <- df[train_ind, ]
test <- df[-train_ind, ]

PCA

     Następnie dane zostaną zredukowane w oparciu o analizę PCA, co pozwoli na ograniczenie liczby zmiennych. Pozwoli to na zmniejszenie czasu potrzebnego na oszacowanie modelu oraz możliwe, że przyczyni się do redukcji błędów szacowanych modeli.

#df1 <- data.matrix(df)  
pca <- princomp(~ tartar + gender + age + height.cm. + 
                  weight.kg. + waist.cm. + eyesight.left. + 
                  eyesight.right. + hearing.left. + hearing.right. + 
                  systolic + relaxation + fasting.blood.sugar + 
                  Cholesterol + triglyceride + HDL + LDL + 
                  hemoglobin + Urine.protein + serum.creatinine + AST +
                  ALT + Gtp + dental.caries, 
                data = trainset, cor = T)

plot((pca$sdev)^2, type = "c")

     Na podstawie powyższego wykresu można stwierdzić, że optymalna liczba nowych zmiennych (komponentów) powinna być określona na poziomie 10. W celu dokłądniejszego określenia liczby zmiennych poniżej zostały zaprezentowane dodatkowa analiza:

summary(pca)

## Importance of components:
##                           Comp.1    Comp.2     Comp.3     Comp.4     Comp.5
## Standard deviation     2.2271892 1.5566911 1.37652751 1.21827695 1.20618216
## Proportion of Variance 0.2066822 0.1009703 0.07895117 0.06184161 0.06061981
## Cumulative Proportion  0.2066822 0.3076525 0.38660363 0.44844524 0.50906505
##                            Comp.6     Comp.7     Comp.8     Comp.9    Comp.10
## Standard deviation     1.13779722 1.08377979 1.07459906 1.01130687 0.94828080
## Proportion of Variance 0.05394094 0.04894078 0.04811513 0.04261423 0.03746819
## Cumulative Proportion  0.56300599 0.61194676 0.66006189 0.70267613 0.74014431
##                           Comp.11    Comp.12    Comp.13   Comp.14    Comp.15
## Standard deviation     0.91278349 0.90673550 0.87073605 0.7973892 0.77764686
## Proportion of Variance 0.03471557 0.03425705 0.03159089 0.0264929 0.02519728
## Cumulative Proportion  0.77485988 0.80911694 0.84070782 0.8672007 0.89239799
##                           Comp.16    Comp.17    Comp.18    Comp.19     Comp.20
## Standard deviation     0.72636816 0.70649209 0.68324994 0.61029718 0.481687383
## Proportion of Variance 0.02198378 0.02079713 0.01945127 0.01551928 0.009667614
## Cumulative Proportion  0.91438177 0.93517890 0.95463017 0.97014945 0.979817063
##                            Comp.21     Comp.22     Comp.23      Comp.24
## Standard deviation     0.456534598 0.428185901 0.296435712 0.0689156225
## Proportion of Variance 0.008684327 0.007639299 0.003661422 0.0001978901
## Cumulative Proportion  0.988501389 0.996140688 0.999802110 1.0000000000

    Liczba zmiennych PCA poprzez analizę cumulative proportion została określona na poziomie 11 i właśnie dla takiej liczby zostały stworzene nowe zbiory test i train:

transformed_test <- as.data.frame(predict(pca, test)[,1:12])

new_train <- as.data.frame(cbind(trainset$smoking, pca$scores[, 1:12]))
colnames(new_train)[1] <- "smoking"

Drzewa decyzyjne

    W tej części zostaną stworzone modele drzew decyzyjnych dla danych pierwotnych i PCA.

Pierwotne dane

    Stworzenie modelu drzewa decyzyjnego dla danych pierwotnych:

drzewo <- tree(as.factor(smoking)~., data = trainset)     

error <- sum(test$smoking != predict(drzewo, test, type = "class"))
relative_error <- (error/nrow(test))
round(1 - relative_error, 4)

## [1] 0.7254

    Błąd dla tego modelu wynosi 0.7279, co wskazuje, że model ten w 72,79% poprawnie przewiduje czy dana osoba jest palaczem na podstawie dostępnych zmiennych.

plot(drzewo)
text(drzewo)

PCA dane

Stworzenie modelu drzewa decyzyjnego dla danych PCA:

drzewo2 <- tree(as.factor(smoking)~., data = new_train)     

error2 <- sum(test$smoking != predict(drzewo2, transformed_test, type = "class"))
relative_error2 <- (error2/nrow(transformed_test))
round(1 - relative_error2, 4)

## [1] 0.6988

    Błąd dla tego modelu wynosi 0.6988, co wskazuje, że model ten w 69,88% poprawnie przewiduje czy dana osoba jest palaczem na podstawie dostępnych zmiennych. Jest on więc gorzej dopasowany niż model dla danych pierwotnych.

plot(drzewo2)
text(drzewo2)

Neural networks

    Na tym etapie zostaną stworzone dwa modele wykorzystujące sieci neuronowe dla danych pierwotnych oraz zreukowanych po PCA. Najpierw natomiast została stworzona funkcja optimal_size, która miała za zadanie sprawdzić dla jakiej ilości neuronów nasz model jest najlepiej dopasowany. W przypadku obu najlepiej dopasowały się dla 6 neuronów.

optimal_size <- function(dane, test1){
  bledy <- list()
  for (i in 2:10){
    blad = 0
    model = nnet(dane, class.ind(smoking), softmax = T, size = i, maxit = 1000, decay = 0.1)
    blad = sum(predict(model, test1, type="class") != test$smoking)
    bledy[i] <- blad
    }
    
  min_value <- min(unlist(bledy))
  min_index <- match(min_value, bledy)
  return(min_index)
}
# optimal_size(dane, test)               wynik optymalizacji, 6 neuronóW
# optimal_size(dane2, transformed_test)  dla danych PCA najlepsza wartość size to również 6

Pierwotne dane

attach(trainset)
dane <- cbind(tartar, gender, age, height.cm.,  
                  weight.kg., waist.cm., eyesight.left., 
                  eyesight.right., hearing.left., hearing.right., 
                  systolic, relaxation, fasting.blood.sugar, 
                  Cholesterol, triglyceride, HDL, LDL, 
                  hemoglobin, Urine.protein, serum.creatinine, AST,
                  ALT, Gtp, dental.caries)

model <- nnet(dane, class.ind(smoking), softmax = T, size = 6, maxit = 1000)
blad <- sum(predict(model, test, type="class")!= test$smoking)

round(1 - (blad/nrow(test)),4)

## [1] 0.7208

Dopasowanie dla danych pierwotnych wynosi ponad 72%

PCA dane

attach(new_train)
dane2 <- cbind(Comp.1, Comp.2, Comp.3, Comp.4, Comp.5, 
              Comp.6, Comp.7, Comp.8, Comp.9, Comp.10, 
              Comp.11, Comp.12)

model2 <- nnet(dane2, class.ind(smoking), softmax=T, size=6, maxit=1000, decay = 0.05)
blad2 <- sum(predict(model2, transformed_test, type="class")!= test$smoking)

round(1-(blad2/nrow(transformed_test)),4)

## [1] 0.7528

    Dopasowanie modelu sieci neuronowych dla danych PCA wynosi ponad 75% co wsazuje, że jest to najlepszy najbardziej trafjny model spośród wszystkich badanych.

Podsumowanie

     Podsumowując, w projekcie zostały zaprezentowane 4 modele tworzone na podstawie drzewa decyzyjnego oraz sieci neuronowych z wykorzystaniem danych zwykłych oraz przetworzonych za pomocą metody PCA. Z analizy błędu można stwierdzić, że najlepiej dopasowanym modelem jest ostatni badany, który wykorzystuje sieci neuronowe oraz dane PCA. Dopasowanie dla niego wynosi ponad 75%. Jesteśmy więc w stanie stwierdzić, że z 72% dokładnością możemy określić na podstawie “sygnałów” ciała czy dana osoba jest palaczem.
     Drzewa decyzyjne z kolei pokazały, że jednym z kluczowych elementów branych pod uwagę jest płeć. Wynika to z tego, że statystycznie więcej mężczyzn jest palaczami, co zwiększa prawdopodobieńśtwo, że właśnie osoba z tej grupy będzie palaczem. Jeżeli chcialibyśmy rzeczywiście działać tylko i wyłącznie na danych dotyczących organizmu powinniśmy pozbyć się zmiennych takich jak: wzrost, płeć, wiek itp., które bardziej określają cechy osób.