Assigment 2: UCI Dataset
1.Import the data to R.
dades <- read.csv("C:/Users/HP/Documents/ICU.csv", header = T)Fem un attach per poder utilitzar les variables del dataframe en el workspace i miro les característiques de les dades
attach(dades)
class(dades)## [1] "data.frame"dim(dades)## [1] 15690 732.Does the dataset contain more than one entry per Subject?
Per saber si hi ha més d’una entrada per subjecte, miro si hi han duplicats de la columna Subjects. Fent table es pot veure el número de duplicats TRUE i FALSE.
table(duplicated(dades[,1]))##
## FALSE
## 15690Ha resultat que no hi ha cap duplicat en la columna de subjects, per tant, només hi ha una entrada per a cada subjecte
3.Convert your variables to its correct format (factors, numeric, etc).
S’ha de tenir en compte que es tracta amb diferents tipus de dades, entre les quals hi ha factors (0 i 1 o 1 i 2). Així doncs, es mira el dataset per veure quines variables són i posar-les a factor. S’ha mirat a la literatura que volen dir cada un dels números abans de passar-los a factor.
Primer miro quina classe és cada variable
vf <- names(dades)
for (i in vf){print (class(dades[,i]))}## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "numeric"
## [1] "integer"
## [1] "integer"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "integer"
## [1] "integer"
## [1] "numeric"
## [1] "integer"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "integer"
## [1] "numeric"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "numeric"
## [1] "numeric"
## [1] "numeric"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "integer"
## [1] "numeric"
## [1] "integer"
## [1] "integer"
## [1] "numeric"
## [1] "integer"
## [1] "numeric"A continuació canvio les que haurien de ser factors
vf <- names(dades)[23:50]
for (f in vf) {dades[,f] <- as.factor(dades[,f])}
for (i in vf) {dades[,i] <- factor(c("neg","pos")[c(dades[,i])])}
vf <- names(dades)[54:58]
for (f in vf) {dades[,f] <- as.factor(dades[,f])}
for (i in vf) {dades[,i] <- factor(c("neg","pos")[c(dades[,i])])}
# D'aquestes columnes no modifico les variables perquè després s'utilitzaran. Nomes les passo a factor
vf <- names(dades)[59:61]
for (f in vf) {dades[,f] <- as.factor(dades[,f])}
#Gènere
dades[,62] <- as.factor(c(dades[,62]))
dades[,62] <- as.factor(c("fem","masc")[c(dades[,62])])
#Obesitat
dades[,64] <- as.factor(c(dades[,64]))
dades[,64] <- as.factor(c("neg","pos")[c(dades[,64])])4. Use an outlier detector and remove odd/abnormal samples.
Primer de tot, per treure els outliers s’instal•la la llibreria “outliers”. La funció de detecció dels outliers no pot treballar amb factors així que se seleccionen només les variables numeric i integer. Tampoc s’agafen les 3 primeres columnes perquè són identificadors
A continuació s’estableix un p-valor a partir del qual escollirà si el punt es considera outlier o no. En aquest cas s’ha escollir p-value < 0.05.
Per comprovar si s’eliminen els outliers es fa un boxplot abans i després de treure’ls
boxplot(dades[15:20], main='Boxplot abans de treure outlier')
Ara s’utilitza la llibreria outlier
library(outliers)
for (i in 4:length(dades)) {
if (class(dades[,i]) == 'numeric' || class(dades[,i]) == 'integer'){
a <- dades[,i]
resul <- chisq.out.test(a, variance=var(a), opposite = FALSE)
resul2 <- chisq.out.test(a, variance=var(a), opposite = TRUE)
if ((resul$p.value<0.05 & is.na(resul$p.value)==FALSE) | (resul2$p.value<0.05) & (is.na(resul2$p.value)==FALSE)) {
b <- (rm.outlier(a, fill = TRUE, median = FALSE, opposite = FALSE))
dades[,i] <- b
}
}
}Mirant els dos boxplot es pot veure com, efectivament, s’han tret valors outliers d’aquestes variables.
boxplot(dades[15:20], main='Boxplot després de treure outliers')
A continuació, perquè l’anàlisi de les dades sigui millor, faig una tipificació. és a dir, faig que totes les variables que no siguin factor tinguin mitjana 0 i SD 1. Tampoc incloc les 3 primeres columnes d’identificadors
for (i in 4:length(dades)) {
if (class(dades[,i]) == 'numeric' || class(dades[,i]) == 'integer'){
dades[,i] <- (dades[,i]-mean(dades[,i]))/sd(dades[,i])
}
}
#Comprovo que s'ha fet bé
for (i in 4:length(dades)) {
if (class(dades[,i]) == 'numeric' || class(dades[,i]) == 'integer'){
print(sd(dades[,i]))
print(mean(dades[,i]))
}
}## [1] 1
## [1] 6.48021e-17
## [1] 1
## [1] -3.654584e-18
## [1] 1
## [1] 4.340609e-17
## [1] 1
## [1] -3.458442e-16
## [1] 1
## [1] 2.393956e-16
## [1] 1
## [1] 1.080552e-16
## [1] 1
## [1] 1.021883e-16
## [1] 1
## [1] 1.0711e-16
## [1] 1
## [1] 2.422881e-16
## [1] 1
## [1] 1.831183e-16
## [1] 1
## [1] -8.539256e-17
## [1] 1
## [1] 2.120975e-16
## [1] 1
## [1] -3.934146e-16
## [1] 1
## [1] 7.230182e-17
## [1] 1
## [1] 5.848573e-17
## [1] 1
## [1] -3.161537e-16
## [1] 1
## [1] 1.705944e-17
## [1] 1
## [1] 1.686226e-16
## [1] 1
## [1] -1.610896e-17
## [1] 1
## [1] 9.935047e-18
## [1] 1
## [1] -4.707144e-17
## [1] 1
## [1] 3.232256e-18
## [1] 1
## [1] 9.346617e-17
## [1] 1
## [1] -4.948607e-16
## [1] 1
## [1] -3.3497e-17
## [1] 1
## [1] 5.778521e-18
## [1] 1
## [1] -4.542696e-17
## [1] 1
## [1] 2.23681e-17
## [1] 1
## [1] -1.978299e-17
## [1] 1
## [1] -3.157285e-17
## [1] 1
## [1] 1.75469e-16
## [1] 1
## [1] 9.655248e-17Fent print de la mitjana i la SD de cada variable no factor es pot observar que totes tenen SD=1 i una mean molt propera a zero (~e-17).
5.Find the variables related with the Simplied Acute Physiology Score (SAPS) and the Sequential Organ Failure Assessment (SOFA). Explain both scores.
SAPS i SOFA són uns paràmetres que ens indiquen la severitat del pacient en la UCI basant-se en diferents paràmetres. Concretament, el SAPS considera els paràmetres següents: * Edat * Freqüència cardíaca (HR) * Pressió Sistòlica * Temperatura * Escala de coma de Glasgow * Ventilació mecànica o CPAP * PaO2 * FiO2 * Orina * Nitrogen en urea (BUN) * Sodiu * Potassi * Bicarbonat * Bilirubina * Glòbuls blancs * Malalties cròniques * Tipus d’admissió
D’altra banda, el SOFA utilitza els paràmetres següents, considerant el sistema cardíac, nerviós, renal i fetge: * PaO2/FiO2 (mmHg) * Escala de coma de Glasgow * Pressió mitjana arterior o administració de vasopressors * Bilirubina (mg/dl)[µmol/L] * Plaquetes x103/µl * Creatinina o orina output (mg/dl) [µmol/L]
6.Build two prediction models, each one predicting the SOFA and SAPS indices.
En aquest apartat es crearan dos models que prediguin els índex SAPS i SOFA. Es faran aquests dos models de dues maneres i després es comprovarà quin és el millor mètode. Per a fer això, es fa una cross-validation de les dades, és a dir, una meitat de les dades serviran per fer el model, i la altra s’utilitzaran en la validació. Cal destacar que hi ha una variable que és sempre 0, així que s’elimina perquè no provoqui problemes en fer el model.
A continuació se separen les dues meitats.
dades2 <- dades[1:(nrow(dades)/2),]
#Treiem BLOOD perquè és un factor que només té un valor
dades2 <- dades2[,-45]
#2a meitat que s'utilitzarà per la validació
dades3=dades[((nrow(dades)/2)+1):nrow(dades),]
dades3 <- dades3[,-45]Així doncs, primer es fa el model de l’índex SAPS. S’utilitza la funció step, que implementa un model de regressió tenint en compte un conjunt de variables. S’especifica un conjunt de variables ampli (anomenat upper) i un conjunt de variables més petit (lower), i ambdós conjunts formaran l’scope de la funció.Aleshores, la funció anirà del conjunt més petit al més gran (direction=forward), del més gran al més petit (backward), o en ambdós sentits. En aquest procés triarà les variables que donin lloc a un model adequat per definir la variable y que se li ha especificat. Així doncs, step intenta trobar la combinació de variables que donen un valor de AIC més òptim.
En el meu cas, he triat una direcció de forward, del conjunt petit al gran.
#Posem el nivell més baix només incloent l'edad
lm_SAPS_lower <- lm(ICUSTAY_ADMIT_SAPS~ICUSTAY_ADMIT_AGE, data = dades2)
#Posem el nivell més alt incloent tots menys identificadors i SAPS i SOFA
lm_SAPS_upper <- lm(ICUSTAY_ADMIT_SAPS~.-SUBJECT_ID-HADM_ID-ICUSTAY_ID-ICUSTAY_ADMIT_SOFA-ICUSTAY_ADMIT_SAPS, data = dades2)
step_SAPS <- step(lm_SAPS_lower,scope=list(lower=lm_SAPS_lower, upper=lm_SAPS_upper), direction='forward')*S’obté un model amb un AIC=-7312.99. Com més petit aquest valor AIC millor.
Els paràmeters utilitzats per crear el model són els següents:
ICUSTAY_ADMIT_SAPS ~ ICUSTAY_ADMIT_AGE + mech_vent_num + Bicarbmin + HR_DELTA + WBC_MAX + Creatmax + WBC_MIN + TEMP_MAX + HR_MIN + HR_MAX + vasopress_num + mincumvol + oneyear_num + gender_num + TEMP_MIN + FLUID_BALANCE + Bcmaxratio + hospital_mort_num + DIABETES_COMPLICATED + WEIGHT_ADMIT + maxcumvol + DIABETES_UNCOMPLICATED + hemo_num + NUM_SIRS_COND + sirs_num + PERIPHERAL_VASCULAR + CARDIAC_ARRHYTHMIAS + CR_ADMIT + METASTATIC_CANCER + ELIX_HOSPITAL_MORT_PT + HYPOTHYROIDISM + WBC_DISCH + avghr + stddevhr + LIVER_DISEASE
Dels paràmetres que han sortit significants, n’hi ha que coincideixen amb els que posteriorment s’utilitzaran en el model fet seguint la literatura, com l’edat, el heart rate, la temperatura, la concentració de bicarbonat, l’orina, i també malalties cròniques com la diabetis, cancer o problemes de fetge, cardíacs i respiratoris. Destacar que amb l’step surt significància amb el gènere o el pes, paràmetres que no són contemplats en els models de SAPS i SOFA.
A continuació es fa el model tenint en compte les variables que especifica el SAPS. Com s’ha comentat, n’hi ha que coincideixen amb les trobades per l’step i n’hi ha que no.
lm_SAPS <- lm(ICUSTAY_ADMIT_SAPS~ ICUSTAY_ADMIT_AGE+HR_ADMIT+HR_DELTA+HR_DISCH+HR_MAX+HR_MIN+TEMP_MAX+TEMP_MIN+URINE_TOTAL+mech_vent_num+Bicarbmin+BUNmax+HYPOTHYROIDISM+PSYCHOSES+CHRONIC_PULMONARY+PARALYSIS+AIDS+CONGESTIVE_HEART_FAILURE+LYMPHOMA+HYPERTENSION+RENAL_FAILURE+LIVER_DISEASE+METASTATIC_CANCER+SOLID_TUMOR, data = dades2)
summary(lm_SAPS)##
## Call:
## lm(formula = ICUSTAY_ADMIT_SAPS ~ ICUSTAY_ADMIT_AGE + HR_ADMIT +
## HR_DELTA + HR_DISCH + HR_MAX + HR_MIN + TEMP_MAX + TEMP_MIN +
## URINE_TOTAL + mech_vent_num + Bicarbmin + BUNmax + HYPOTHYROIDISM +
## PSYCHOSES + CHRONIC_PULMONARY + PARALYSIS + AIDS + CONGESTIVE_HEART_FAILURE +
## LYMPHOMA + HYPERTENSION + RENAL_FAILURE + LIVER_DISEASE +
## METASTATIC_CANCER + SOLID_TUMOR, data = dades2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.52314 -0.43186 -0.01807 0.41895 3.13704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.618095 0.015788 -39.151 < 2e-16 ***
## ICUSTAY_ADMIT_AGE 0.290266 0.008573 33.857 < 2e-16 ***
## HR_ADMIT 0.009203 0.011659 0.789 0.429898
## HR_DELTA 0.041335 0.075102 0.550 0.582072
## HR_DISCH -0.012689 0.061610 -0.206 0.836833
## HR_MAX 0.197276 0.058086 3.396 0.000686 ***
## HR_MIN -0.163573 0.010636 -15.380 < 2e-16 ***
## TEMP_MAX 0.119534 0.008268 14.457 < 2e-16 ***
## TEMP_MIN -0.054601 0.008174 -6.680 2.56e-11 ***
## URINE_TOTAL -0.026442 0.007939 -3.331 0.000870 ***
## mech_vent_numpos 0.976936 0.017041 57.328 < 2e-16 ***
## Bicarbmin -0.170242 0.008533 -19.952 < 2e-16 ***
## BUNmax 0.118194 0.009088 13.006 < 2e-16 ***
## HYPOTHYROIDISMpos 0.099740 0.029388 3.394 0.000692 ***
## PSYCHOSESpos -0.014942 0.043037 -0.347 0.728463
## CHRONIC_PULMONARYpos 0.004429 0.021242 0.209 0.834839
## PARALYSISpos 0.084347 0.067834 1.243 0.213750
## AIDSpos 0.118555 0.120475 0.984 0.325111
## CONGESTIVE_HEART_FAILUREpos -0.014302 0.020421 -0.700 0.483725
## LYMPHOMApos 0.029419 0.058869 0.500 0.617274
## HYPERTENSIONpos -0.016319 0.017008 -0.959 0.337337
## RENAL_FAILUREpos 0.140423 0.036112 3.889 0.000102 ***
## LIVER_DISEASEpos 0.095778 0.038051 2.517 0.011852 *
## METASTATIC_CANCERpos -0.021933 0.037691 -0.582 0.560651
## SOLID_TUMORpos 0.021711 0.025435 0.854 0.393353
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6654 on 7820 degrees of freedom
## Multiple R-squared: 0.5599, Adjusted R-squared: 0.5586
## F-statistic: 414.6 on 24 and 7820 DF, p-value: < 2.2e-16Repetim el procés amb el SOFA
#Posem el nivell més baix només incloent l'edad
lm_SOFA_lower <- lm(ICUSTAY_ADMIT_SOFA~ICUSTAY_ADMIT_AGE, data = dades2)
#Posem el nivell més alt incloent tots menys identificadors i SAPS i SOFA
lm_SOFA_upper <- lm(ICUSTAY_ADMIT_SOFA~.-SUBJECT_ID-HADM_ID-ICUSTAY_ID-ICUSTAY_ADMIT_SOFA-ICUSTAY_ADMIT_SAPS, data = dades2)
step_SOFA <- step(lm_SOFA_lower,scope=list(lower=lm_SOFA_lower, upper=lm_SOFA_upper), direction='forward')*S’obté un model amb AIC=-8662.61
El model creat per l’step té els components següents:
ICUSTAY_ADMIT_SOFA ~ ICUSTAY_ADMIT_AGE + mech_vent_num + vasopress_num + BUNmax + mincumvol + Creatmax + HR_DELTA + COAGULOPATHY + LIVER_DISEASE + NUM_SIRS_COND + Bicarbmin + FLUID_BALANCE + maxcumvol + twenteight_num + URINE_TOTAL + sepsis_num + WBC_MIN + WBC_MAX + HYPERTENSION + TEMP_MAX + WEIGHT_ADMIT + CR_ADMIT + stddevhr + numoverload + HR_MAX + HR_MIN + AIDS + TEMP_MIN + PERIPHERAL_VASCULAR + CARDIAC_ARRHYTHMIAS + CONGESTIVE_HEART_FAILURE + HR_ADMIT + Bcmaxratio + hemo_num + gender_num + OTHER_NEUROLOGICAL + RENAL_FAILURE + OBESITY + oneyear_num
Com en el cas del SAPS, dels paràmetres que han sortit significants, n’hi ha que coincideixen amb els que posteriorment s’utilitzaran en el model fet seguint la literatura, com l’edat, el heart rate, la temperatura, la concentració de bicarbonat, l’orina, i també malalties cròniques com la diabetis, cancer o problemes de fetge, cardíacs i respiratoris. Tanmateix, amb el SOFA, es contempla també el paràmetre d’obesitat, juntament amb el gènere i l’edat. Comentar això perquè posteriorment s’estudiarà la relació de l’obsesitat amb un paràmeter que indica la mort o la supervivència del pacient.
A continuació es fa el model tenint en compte les variables que especifica el SAPS. Com s’ha comentat, n’hi ha que coincideixen amb les trobades per l’step i n’hi ha que no.
lm_SOFA <- lm(ICUSTAY_ADMIT_SOFA~ mech_vent_num+PULMONARY_CIRCULATION+vasopress_num+HYPERTENSION+LIVER_DISEASE+COAGULOPATHY+DEFICIENCY_ANEMIAS+BUNmax+URINE_TOTAL+RENAL_FAILURE, data = dades2)
summary(lm_SOFA)##
## Call:
## lm(formula = ICUSTAY_ADMIT_SOFA ~ mech_vent_num + PULMONARY_CIRCULATION +
## vasopress_num + HYPERTENSION + LIVER_DISEASE + COAGULOPATHY +
## DEFICIENCY_ANEMIAS + BUNmax + URINE_TOTAL + RENAL_FAILURE,
## data = dades2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9951 -0.3244 -0.0387 0.3553 3.1794
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.933411 0.013976 -66.788 < 2e-16 ***
## mech_vent_numpos 1.035908 0.016460 62.935 < 2e-16 ***
## PULMONARY_CIRCULATIONpos 0.008961 0.053414 0.168 0.86678
## vasopress_numpos 0.626397 0.016166 38.748 < 2e-16 ***
## HYPERTENSIONpos -0.083378 0.015717 -5.305 1.16e-07 ***
## LIVER_DISEASEpos 0.473924 0.035722 13.267 < 2e-16 ***
## COAGULOPATHYpos 0.362578 0.029328 12.363 < 2e-16 ***
## DEFICIENCY_ANEMIASpos -0.003739 0.023111 -0.162 0.87146
## BUNmax 0.253635 0.007715 32.876 < 2e-16 ***
## URINE_TOTAL -0.051218 0.007322 -6.995 2.87e-12 ***
## RENAL_FAILUREpos 0.099711 0.033976 2.935 0.00335 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6299 on 7834 degrees of freedom
## Multiple R-squared: 0.6094, Adjusted R-squared: 0.6089
## F-statistic: 1222 on 10 and 7834 DF, p-value: < 2.2e-16Es poden mirar es coeficients dels quatre models
coef(step_SAPS)## (Intercept) ICUSTAY_ADMIT_AGE
## -0.54840051 0.27389546
## mech_vent_numpos Bicarbmin
## 0.87152677 -0.09983670
## HR_DELTA WBC_MAX
## 0.03974214 0.29782412
## Creatmax WBC_MIN
## 0.08746188 -0.19219288
## TEMP_MAX HR_MIN
## 0.09609597 -0.18326649
## HR_MAX vasopress_numpos
## 0.17510665 0.15035082
## mincumvol oneyear_num2
## 0.10300879 0.09533236
## gender_nummasc TEMP_MIN
## -0.07187130 -0.04480508
## FLUID_BALANCE Bcmaxratio
## -0.11766048 0.03185776
## hospital_mort_num2 DIABETES_COMPLICATEDpos
## 0.11894997 0.13146457
## WEIGHT_ADMIT maxcumvol
## -0.02944028 0.07303371
## DIABETES_UNCOMPLICATEDpos hemo_numpos
## 0.05559112 0.12471844
## NUM_SIRS_COND sirs_numpos
## 0.03999319 -0.07827330
## PERIPHERAL_VASCULARpos CARDIAC_ARRHYTHMIASpos
## 0.06502252 -0.05640668
## CR_ADMIT METASTATIC_CANCERpos
## 0.03600222 -0.10335000
## ELIX_HOSPITAL_MORT_PT HYPOTHYROIDISMpos
## 0.01962061 0.05386466
## WBC_DISCH avghr
## 0.01666831 0.02622099
## stddevhr LIVER_DISEASEpos
## -0.01621427 0.05509361coef(lm_SAPS)## (Intercept) ICUSTAY_ADMIT_AGE
## -0.618095371 0.290265896
## HR_ADMIT HR_DELTA
## 0.009203303 0.041334643
## HR_DISCH HR_MAX
## -0.012688749 0.197276427
## HR_MIN TEMP_MAX
## -0.163573140 0.119533742
## TEMP_MIN URINE_TOTAL
## -0.054600699 -0.026441582
## mech_vent_numpos Bicarbmin
## 0.976936471 -0.170242283
## BUNmax HYPOTHYROIDISMpos
## 0.118194387 0.099740330
## PSYCHOSESpos CHRONIC_PULMONARYpos
## -0.014941734 0.004429007
## PARALYSISpos AIDSpos
## 0.084346677 0.118555483
## CONGESTIVE_HEART_FAILUREpos LYMPHOMApos
## -0.014301660 0.029419061
## HYPERTENSIONpos RENAL_FAILUREpos
## -0.016319051 0.140423081
## LIVER_DISEASEpos METASTATIC_CANCERpos
## 0.095778247 -0.021932548
## SOLID_TUMORpos
## 0.021711026coef(step_SOFA)## (Intercept) ICUSTAY_ADMIT_AGE
## -0.89484282 0.05053986
## mech_vent_numpos vasopress_numpos
## 0.94019760 0.52146574
## BUNmax mincumvol
## 0.12879146 0.14037056
## Creatmax HR_DELTA
## 0.02569621 0.04649421
## COAGULOPATHYpos LIVER_DISEASEpos
## 0.28991462 0.36629269
## NUM_SIRS_COND Bicarbmin
## 0.04474599 -0.05914110
## FLUID_BALANCE maxcumvol
## -0.26313273 0.24262376
## twenteight_num2 URINE_TOTAL
## 0.14487110 -0.05590902
## sepsis_numpos WBC_MIN
## 0.13989073 -0.16094628
## WBC_MAX HYPERTENSIONpos
## 0.16254951 -0.08017517
## TEMP_MAX WEIGHT_ADMIT
## 0.04413850 0.02113932
## CR_ADMIT stddevhr
## 0.07831984 -0.04490322
## numoverload HR_MAX
## -0.02804245 0.04648790
## HR_MIN AIDSpos
## -0.04406294 0.27528754
## TEMP_MIN PERIPHERAL_VASCULARpos
## -0.01937373 0.06069618
## CARDIAC_ARRHYTHMIASpos CONGESTIVE_HEART_FAILUREpos
## -0.05664523 0.05257634
## HR_ADMIT Bcmaxratio
## 0.02232735 -0.02172209
## hemo_numpos gender_nummasc
## 0.10360891 0.03086253
## OTHER_NEUROLOGICALpos RENAL_FAILUREpos
## 0.06890090 -0.05816997
## OBESITYpos oneyear_num2
## 0.02808637 0.03321171coef(lm_SOFA)## (Intercept) mech_vent_numpos PULMONARY_CIRCULATIONpos
## -0.933411477 1.035907522 0.008960778
## vasopress_numpos HYPERTENSIONpos LIVER_DISEASEpos
## 0.626397125 -0.083378429 0.473923881
## COAGULOPATHYpos DEFICIENCY_ANEMIASpos BUNmax
## 0.362578337 -0.003739484 0.253635370
## URINE_TOTAL RENAL_FAILUREpos
## -0.051218491 0.099710643Cal destacar que a l’hora de fer els models basant-se en la literatura, s’han escollit molts més paràmetres pel model SAPS que pel model SOFA (25 VS 11). Això pot tenir un efecte pel que fa al poder predictiu del model. D’altra banda, la funció step ha trobat aproximadament el mateix número de paràmetres pel SAPS i pel SOFA, 36 i 40, respectivament.
Finalment es pot comprovar quin dels dos models és millor, el de l’step o el que es basa en les dades de la literatura. S’utilitza la funció predict i la meitat de les dades que no s’han utilitzat per fer els models
predict_SAPS <- predict(step_SAPS, dades3)
summary(predict_SAPS)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.88800 -0.53800 0.09580 0.03843 0.56470 6.27200predict_SAPS2 <- predict(lm_SAPS, dades3)
summary(predict_SAPS2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.85500 -0.52910 0.10680 0.02539 0.53330 2.55900#Repetim el procés per SOFA
predict_SOFA <- predict(step_SOFA, dades3)
summary(predict_SOFA)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.60400 -0.69070 0.11680 0.04095 0.60680 3.93300predict_SOFA2 <- predict(lm_SOFA, dades3)
summary(predict_SOFA2)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.35900 -0.70160 0.05068 0.03445 0.61230 2.88900#Fem els plots comparant els quatre models
par(mfrow=c(1,2))
plot(predict_SAPS, ylab='prediction y step', xlab='number of patients', col='blue')
plot(predict_SAPS2, ylab='prediction y SAPS variables', xlab='number of patients', col='red')
par(mfrow=c(1,2))
plot(predict_SOFA, ylab='prediction y step', col='blue')
plot(predict_SOFA2, ylab='prediction y SOFA variables', col='red')
dev.off()## null device
## 1Podem observar que per ambdos casos el model que utilitza les variables establertes dóna resultats més aproximats, ja que es veu una dispersió menor. A part d’això, la predicció del índex SOFA per step és millor que la del índex SAPS, mentre que sembla que per la predicció amb les variables, el model SAPS és més acurat. S’ha de tenir en compte però, que al fer els models amb les variables, s’han escollit les variables de manera subjectiva, pensant quines estaven relacionades amb els paràmetres que indicaven en la literatura, però potser en falten.
Tanmateix, mirant els valors del Residual standard error o els Multiple i Adjusted R-square, podem concloure que, en general, no són models molt ajustats.
7.Create a unique numerical variable representing the survival from the oneyear_num (1 year survival), twenteight_num (28 days survival) and hospital_mort_num (in hostpital survival).
En les dades hi ha 3 variables que indiquen si el pacient ha mort o no: oneyear_num (1 year survival), twenteight_num (28 days survival) and hospital_mort_num (in hostpital survival). Si el pacient mor es posa un 2, i si sobreviu, un 1. La idea es crear una altra variable que inclogui la informació de la supervivència. Així doncs, es sumen aquestes 3 columnes i es crea una de nova (indexmort), que tindrà: * 3: no sobreviu * 2: sobrevius fins 28 dies * 1: sobreviu fins 1 any * 0: sobreviu
#Passo les variables de factor a numèric per poder sumar les columnes
vf <- names(dades)[59:61]
for (i in vf) {dades[,i] <- as.numeric(c(1,2)[c(dades[,i])])}
indexmort <- ((dades[,59] + dades[,60] + dades[,61])-3)
dades <- as.vector(cbind(dades, indexmort))
#Miro quans pacients hi ha de cada
table(dades[,74])##
## 0 1 2 3
## 11814 1587 515 17748.Test the relationship between obesity and 1 year survival. Interpret the results.
Per mirar si hi ha relació entre ambdues variables fa una taula de contingència i es mira si surt un p-value significatiu.
table(dades$OBESITY, dades$oneyear_num)##
## 1 2
## neg 7691 2807
## pos 4123 1069summary(table(OBESITY, oneyear_num))## Number of cases in table: 15690
## Number of factors: 2
## Test for independence of all factors:
## Chisq = 70.62, df = 1, p-value = 4.341e-17En aquest cas ha donat un p-value significatiu (p-value = 1.149e-08), per tant, són dues variables que tenen una dependència. Considerant que negatiu és no tenir Obesitat i 1 indica que sobreviu, es pot veure que sobreviuen més les persones amb obesitat que les que no. És al contrari que podríem esperar reflexionant sobre la relació entre aquestes dues variables.
9.Split randomly your data in two datasets. Build a prediction model for explaining the survival with the first dataset.
Ara tornem a partir les dades en dues meitats per poder fer la cross-validation del model. Ara també inclouran la columna de indexmort. Abans de tot, converteixo a factor un altre cop les variables que s’han passat a numèric en l’anterior apartat 7.
vf <- names(dades)[59:61]
for (f in vf) {dades[,f] <- as.factor(dades[,f])}
#Ara separem en dos meitat
dades2 <- dades[1:(nrow(dades)/2),]
#Treiem BLOOD perquè és un factor que només té un valor
dades2 <- dades2[,-45]
#2a meitat que s'utilitzarà per la validació
dades3=dades[((nrow(dades)/2)+1):nrow(dades),]
dades3 <- dades3[,-45]A continuació fem servir step per fer el model de predicció de supervivència, d’una manera semblant com s’ha fet els de SAPS i SOFA.
#Posem el nivell més baix només incloent l'edad
lm_surv_lower <- lm(indexmort~ICUSTAY_ADMIT_AGE, data = dades2)
#Posem el nivell més alt incloent tots menys identificadors, SAPS i SOFA i mortalitats i algunes variables que es creu que no tenen relació
lm_surv_upper <- lm(indexmort~.-SUBJECT_ID-HADM_ID-ICUSTAY_ID-ICUSTAY_ADMIT_SOFA-ICUSTAY_ADMIT_SAPS-indexmort-oneyear_num-hospital_mort_num-twenteight_num, data = dades2)
step_surv <- step(lm_surv_lower,scope=list(lower=lm_surv_lower, upper=lm_surv_upper), direction='forward')*Aquest model té un AIC=-3599.97. El valor és més gran i, per tant, pitjor, que els del models de SAPS i SOFA.
#Els coeficients del step són els següents:
coef(step_surv)## (Intercept) ICUSTAY_ADMIT_AGE
## 0.40107429 0.12255787
## HR_DELTA BUNmax
## -0.02065769 0.18127290
## ELIX_28D_MORT_PT WBC_DISCH
## 0.07829749 0.12328003
## ICUstay_num HR_DISCH
## 0.11458776 -0.31235722
## avghr stddevhr
## 0.25282918 0.09419605
## METASTATIC_CANCERpos sepsis_numpos
## 0.48758413 0.26083710
## URINE_TOTAL WEIGHT_ADMIT
## -0.09993324 -0.04987759
## FLUID_BALANCE mincumvol
## 0.21895860 -0.10766419
## CR_ADMIT HR_ADMIT
## -0.04138544 -0.05293642
## OTHER_NEUROLOGICALpos PULMONARY_CIRCULATIONpos
## 0.23721081 0.23993060
## maxcumvol WBC_MIN
## -0.11942905 -0.04512156
## vasopress_numpos LYMPHOMApos
## 0.07115491 0.25242038
## CONGESTIVE_HEART_FAILUREpos HYPERTENSIONpos
## 0.08124740 -0.04290425
## Bcmaxratio DRUG_ABUSEpos
## 0.03172535 0.12059678
## Bicarbmin AIDSpos
## -0.02261512 0.25796695
## TEMP_MAX DIABETES_UNCOMPLICATEDpos
## -0.01937223 -0.03683655
## mech_vent_numpos HR_MIN
## 0.03286288 -0.02344569Tenint en compte els resultats, step ha utilitzat les variables següents per a fer el model:
indexmort ~ ICUSTAY_ADMIT_AGE + HR_DELTA + BUNmax + ELIX_28D_MORT_PT + WBC_DISCH + ICUstay_num + HR_DISCH + avghr + stddevhr + METASTATIC_CANCER + sepsis_num + URINE_TOTAL + WEIGHT_ADMIT + FLUID_BALANCE + mincumvol + CR_ADMIT + HR_ADMIT + OTHER_NEUROLOGICAL + PULMONARY_CIRCULATION + maxcumvol + WBC_MIN + vasopress_num + LYMPHOMA + CONGESTIVE_HEART_FAILURE + HYPERTENSION + Bcmaxratio + DRUG_ABUSE + Bicarbmin + AIDS + TEMP_MAX + DIABETES_UNCOMPLICATED + mech_vent_num + HR_MIN
Variables com l’edat, el heart rate, la temperatura o la concentració de bicarbonat, ja eren paràmetres que s’esperava que tinguéssin una relació pel que fa a l’índex creat. Tanmateix, cal destacar també que hi ha un gran nombre de malalties amb les quals també hi ha significància. Aquestes malalties acostumen a ser cròniques i, per tant, tindran una implicació important en la supervivència o no dels pacients, per exemple problemes neurològics, cardíacs o respiratòris, diferents tipus de càncer, malalties com la diabetis o l’AIDS, entre d’altres. Destacar també que el fet de consumir drogues també és un paràmetre a tenir en compte.
Així doncs, els paràmetres utilitzats per la funció step es poden entendre des d’un punt de vista mèdic.
Per comprovar el model, es pot fer un predict amb la meitat de les dades no utilitzades per fer el model.
predict_surv <- predict(step_surv, dades3)
summary(predict_surv)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.4790 0.1253 0.3735 0.5153 0.7385 4.3350plot(predict_surv, ylab='prediction y step', main= 'Survival prediction model', col='blue')
dev.off()## null device
## 1Es pot observar que la majoria de punts estan entre -0.5 i 1.5, hi ha una dispersió entre 1.5 i 3, quedant només alguns punts entre 3 i 4.
10.Use the pROC package to plot the receiver operating characteristic (ROC) on your prediction for the 1 year survival and the 28 days survival.
En aquest apartat s’evalua la corba ROC de la predicció feta en funció de dues mesures de supervivència (1 year i 28 days survival).
Que sobrevisquin 1 any significa que indexmort= 1 i sobreviure 28 dies és 2. Per tant, s’escullen els pacients que compleixen aquestes condicions a l’hora de fer la corba ROC.
library(pROC)## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
##
## The following objects are masked from 'package:stats':
##
## cov, smooth, varsurv_one <- dades3$indexmort == 1
proj<-predict(step_surv,dades)
roc(indexmort[surv_one], proj[surv_one], plot=TRUE, main='Corba ROC 1 year survival') 
##
## Call:
## roc.default(response = indexmort[surv_one], predictor = proj[surv_one], plot = TRUE, main = "Corba ROC 1 year survival")
##
## Data: proj[surv_one] in 550 controls (indexmort[surv_one] 0) < 824 cases (indexmort[surv_one] 1).
## Area under the curve: 0.7541#El procés es repeteix per els 28 day survival
surv_28 <- dades3$indexmort == 2
proj<-predict(step_surv,dades)
roc(indexmort[surv_28], proj[surv_28], plot=TRUE, main='Corba ROC 28 days survival') 
##
## Call:
## roc.default(response = indexmort[surv_28], predictor = proj[surv_28], plot = TRUE, main = "Corba ROC 28 days survival")
##
## Data: proj[surv_28] in 182 controls (indexmort[surv_28] 0) < 22 cases (indexmort[surv_28] 1).
## Area under the curve: 0.8482Les corbes ROC observades tenen una AUC de 0.75
library(pls)##
## Attaching package: 'pls'
##
## The following object is masked from 'package:outliers':
##
## scores
##
## The following object is masked from 'package:stats':
##
## loadingsscoreplot(pca, col = as.numeric(dades$indexmort)+2,type = "p", main='PCA Scoreplot')
Fent l’scoreplot de la PCA i separant per colors els índex de mortalitat creats per nosaltres, es pot observar que estan agrupats segons el color. La Cumulative Proportion de la PC1 és petita, al voltant del 9.6%, però té sentit considerant que tenim voltes variables.