Problemas 3.1.
David López
2024-04-18
Caso práctico: Introducción al machine learning con R
#Caso práctico 1. Modelos supervisados. Regresión lineal y múltiple A partir del conjunto de datos Ovarian del paquete survival, haced un estudio de regresión lineal sobre un par de variables que, a priori, penséis que pueden presentar correlación, y un estudio de regresión múltiple seleccionando los mejores predictores.
#Cargamos los datos
library(survival)
data(ovarian)## Warning in data(ovarian): data set 'ovarian' not foundhead(ovarian)## futime fustat age resid.ds rx ecog.ps
## 1 59 1 72.3315 2 1 1
## 2 115 1 74.4932 2 1 1
## 3 156 1 66.4658 2 1 2
## 4 421 0 53.3644 2 2 1
## 5 431 1 50.3397 2 1 1
## 6 448 0 56.4301 1 1 2summary(ovarian)## futime fustat age resid.ds
## Min. : 59.0 Min. :0.0000 Min. :38.89 Min. :1.000
## 1st Qu.: 368.0 1st Qu.:0.0000 1st Qu.:50.17 1st Qu.:1.000
## Median : 476.0 Median :0.0000 Median :56.85 Median :2.000
## Mean : 599.5 Mean :0.4615 Mean :56.17 Mean :1.577
## 3rd Qu.: 794.8 3rd Qu.:1.0000 3rd Qu.:62.38 3rd Qu.:2.000
## Max. :1227.0 Max. :1.0000 Max. :74.50 Max. :2.000
## rx ecog.ps
## Min. :1.0 Min. :1.000
## 1st Qu.:1.0 1st Qu.:1.000
## Median :1.5 Median :1.000
## Mean :1.5 Mean :1.462
## 3rd Qu.:2.0 3rd Qu.:2.000
## Max. :2.0 Max. :2.000#Vemos cuántas variables tienen los datos y cómo son
str(ovarian)## 'data.frame': 26 obs. of 6 variables:
## $ futime : num 59 115 156 421 431 448 464 475 477 563 ...
## $ fustat : num 1 1 1 0 1 0 1 1 0 1 ...
## $ age : num 72.3 74.5 66.5 53.4 50.3 ...
## $ resid.ds: num 2 2 2 2 2 1 2 2 2 1 ...
## $ rx : num 1 1 1 2 1 1 2 2 1 2 ...
## $ ecog.ps : num 1 1 2 1 1 2 2 2 1 2 ...#Observamos correlación entre variables
cor(x=ovarian, method="pearson")## futime fustat age resid.ds rx ecog.ps
## futime 1.00000000 -0.6898795 -0.64836121 -0.38518677 0.24687108 0.01425371
## fustat -0.68987946 1.0000000 0.49586611 0.32433749 -0.15430335 0.22619048
## age -0.64836121 0.4958661 1.00000000 0.29339868 0.04372768 0.12927600
## resid.ds -0.38518677 0.3243375 0.29339868 1.00000000 -0.07784989 -0.14414999
## rx 0.24687108 -0.1543033 0.04372768 -0.07784989 1.00000000 0.00000000
## ecog.ps 0.01425371 0.2261905 0.12927600 -0.14414999 0.00000000 1.00000000#Escogemos las variables futime y age y hacemos el modelo de regresión lineal
regmod<-lm(age~futime, data=ovarian)
summary(regmod)##
## Call:
## lm(formula = age ~ futime, data = ovarian)
##
## Residuals:
## Min 1Q Median 3Q Max
## -18.2440 -3.1173 0.5101 5.2789 15.5201
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 67.723641 3.169255 21.369 < 2e-16 ***
## futime -0.019278 0.004621 -4.172 0.000341 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.848 on 24 degrees of freedom
## Multiple R-squared: 0.4204, Adjusted R-squared: 0.3962
## F-statistic: 17.41 on 1 and 24 DF, p-value: 0.0003408lm(formula = age ~ futime, data = ovarian)##
## Call:
## lm(formula = age ~ futime, data = ovarian)
##
## Coefficients:
## (Intercept) futime
## 67.72364 -0.01928plot(ovarian$futime,ovarian$age, xlab="futime", ylab="age")
abline(regmod)
#Ahora haremos el modelo de regresión lineal múltiple
mrm_ovarian<-lm(ovarian$futime~ovarian$fustat+ovarian$age+ovarian$resid.ds+ovarian$rx+ovarian$ecog.ps, data=ovarian)
summary(mrm_ovarian)##
## Call:
## lm(formula = ovarian$futime ~ ovarian$fustat + ovarian$age +
## ovarian$resid.ds + ovarian$rx + ovarian$ecog.ps, data = ovarian)
##
## Residuals:
## Min 1Q Median 3Q Max
## -374.48 -103.23 -12.18 77.76 384.41
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1276.931 327.157 3.903 0.000882 ***
## ovarian$fustat -309.207 105.682 -2.926 0.008357 **
## ovarian$age -14.373 5.083 -2.828 0.010401 *
## ovarian$resid.ds -48.136 95.331 -0.505 0.619126
## ovarian$rx 125.638 87.254 1.440 0.165364
## ovarian$ecog.ps 109.503 90.484 1.210 0.240315
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 217.3 on 20 degrees of freedom
## Multiple R-squared: 0.6727, Adjusted R-squared: 0.5909
## F-statistic: 8.223 on 5 and 20 DF, p-value: 0.0002355#Seleccionamos los mejores predictores posibles
step(object=mrm_ovarian, direction="both",trace=1)## Start: AIC=285
## ovarian$futime ~ ovarian$fustat + ovarian$age + ovarian$resid.ds +
## ovarian$rx + ovarian$ecog.ps
##
## Df Sum of Sq RSS AIC
## - ovarian$resid.ds 1 12034 956070 283.32
## - ovarian$ecog.ps 1 69130 1013166 284.83
## <none> 944036 285.00
## - ovarian$rx 1 97865 1041901 285.56
## - ovarian$age 1 377400 1321436 291.74
## - ovarian$fustat 1 404069 1348105 292.26
##
## Step: AIC=283.32
## ovarian$futime ~ ovarian$fustat + ovarian$age + ovarian$rx +
## ovarian$ecog.ps
##
## Df Sum of Sq RSS AIC
## <none> 956070 283.32
## - ovarian$ecog.ps 1 88920 1044991 283.64
## - ovarian$rx 1 101328 1057398 283.94
## + ovarian$resid.ds 1 12034 944036 285.00
## - ovarian$age 1 414281 1370351 290.69
## - ovarian$fustat 1 467936 1424007 291.68##
## Call:
## lm(formula = ovarian$futime ~ ovarian$fustat + ovarian$age +
## ovarian$rx + ovarian$ecog.ps, data = ovarian)
##
## Coefficients:
## (Intercept) ovarian$fustat ovarian$age ovarian$rx
## 1213.26 -322.41 -14.82 127.70
## ovarian$ecog.ps
## 120.52#Ejecutamos de nuevo el modelo con los predictores más significativos
mrm_ovarian_predictores<-(lm(ovarian$futime~ovarian$fustat+ovarian$age+ovarian$rx+ovarian$ecog.ps, data = ovarian))
summary(mrm_ovarian_predictores)##
## Call:
## lm(formula = ovarian$futime ~ ovarian$fustat + ovarian$age +
## ovarian$rx + ovarian$ecog.ps, data = ovarian)
##
## Residuals:
## Min 1Q Median 3Q Max
## -377.09 -94.77 -16.40 73.50 400.66
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1213.263 296.477 4.092 0.000521 ***
## ovarian$fustat -322.407 100.565 -3.206 0.004244 **
## ovarian$age -14.824 4.914 -3.017 0.006568 **
## ovarian$rx 127.701 85.599 1.492 0.150606
## ovarian$ecog.ps 120.524 86.240 1.398 0.176842
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 213.4 on 21 degrees of freedom
## Multiple R-squared: 0.6686, Adjusted R-squared: 0.6054
## F-statistic: 10.59 on 4 and 21 DF, p-value: 7.383e-05#Caso práctico 2. Modelos aprendizaje automático y ANOVA A partir del conjunto de datos lung del paquete survival, se pide realizar un breve estudio aplicando los conceptos estudiados en este LAB. Para ello, se pide realizar un estudio de regresión lineal, regresión múltiple, test ANOVA y test clustering sobre dicho conjunto de datos. Las variables elegidas en cada uno de los casos posibles que tratar se dan a elegir al estudiante.
#Cargamos los datos
library(survival)
data(lung)## Warning in data(lung): data set 'lung' not foundhead(lung)## inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1 3 306 2 74 1 1 90 100 1175 NA
## 2 3 455 2 68 1 0 90 90 1225 15
## 3 3 1010 1 56 1 0 90 90 NA 15
## 4 5 210 2 57 1 1 90 60 1150 11
## 5 1 883 2 60 1 0 100 90 NA 0
## 6 12 1022 1 74 1 1 50 80 513 0#Vemos cuántas variables tienen los datos y cómo son
str(lung)## 'data.frame': 228 obs. of 10 variables:
## $ inst : num 3 3 3 5 1 12 7 11 1 7 ...
## $ time : num 306 455 1010 210 883 ...
## $ status : num 2 2 1 2 2 1 2 2 2 2 ...
## $ age : num 74 68 56 57 60 74 68 71 53 61 ...
## $ sex : num 1 1 1 1 1 1 2 2 1 1 ...
## $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
## $ meal.cal : num 1175 1225 NA 1150 NA ...
## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...#Eliminamos los valores NA y hacemos resumen
data_lung=na.omit(lung)
summary(data_lung)## inst time status age
## Min. : 1.00 Min. : 5.0 Min. :1.000 Min. :39.00
## 1st Qu.: 3.00 1st Qu.: 174.5 1st Qu.:1.000 1st Qu.:57.00
## Median :11.00 Median : 268.0 Median :2.000 Median :64.00
## Mean :10.71 Mean : 309.9 Mean :1.719 Mean :62.57
## 3rd Qu.:15.00 3rd Qu.: 419.5 3rd Qu.:2.000 3rd Qu.:70.00
## Max. :32.00 Max. :1022.0 Max. :2.000 Max. :82.00
## sex ph.ecog ph.karno pat.karno
## Min. :1.000 Min. :0.0000 Min. : 50.00 Min. : 30.00
## 1st Qu.:1.000 1st Qu.:0.0000 1st Qu.: 70.00 1st Qu.: 70.00
## Median :1.000 Median :1.0000 Median : 80.00 Median : 80.00
## Mean :1.383 Mean :0.9581 Mean : 82.04 Mean : 79.58
## 3rd Qu.:2.000 3rd Qu.:1.0000 3rd Qu.: 90.00 3rd Qu.: 90.00
## Max. :2.000 Max. :3.0000 Max. :100.00 Max. :100.00
## meal.cal wt.loss
## Min. : 96.0 Min. :-24.000
## 1st Qu.: 619.0 1st Qu.: 0.000
## Median : 975.0 Median : 7.000
## Mean : 929.1 Mean : 9.719
## 3rd Qu.:1162.5 3rd Qu.: 15.000
## Max. :2600.0 Max. : 68.000str(data_lung)## 'data.frame': 167 obs. of 10 variables:
## $ inst : num 3 5 12 7 11 1 7 6 12 22 ...
## $ time : num 455 210 1022 310 361 ...
## $ status : num 2 2 1 2 2 2 2 2 2 2 ...
## $ age : num 68 57 74 68 71 53 61 57 57 70 ...
## $ sex : num 1 1 1 2 2 1 1 1 1 1 ...
## $ ph.ecog : num 0 1 1 2 2 1 2 1 1 1 ...
## $ ph.karno : num 90 90 50 70 60 70 70 80 80 90 ...
## $ pat.karno: num 90 60 80 60 80 80 70 80 70 100 ...
## $ meal.cal : num 1225 1150 513 384 538 ...
## $ wt.loss : num 15 11 0 10 1 16 34 27 60 -5 ...
## - attr(*, "na.action")= 'omit' Named int [1:61] 1 3 5 12 13 14 16 20 23 25 ...
## ..- attr(*, "names")= chr [1:61] "1" "3" "5" "12" ...#Observamos correlación entre variables
cor(x=data_lung, method="pearson")## inst time status age sex
## inst 1.00000000 0.02462876 -0.12719801 0.04859452 0.084362910
## time 0.02462876 1.00000000 -0.16217237 -0.07854153 0.114149616
## status -0.12719801 -0.16217237 1.00000000 0.15911933 -0.218780030
## age 0.04859452 -0.07854153 0.15911933 1.00000000 -0.125280356
## sex 0.08436291 0.11414962 -0.21878003 -0.12528036 1.000000000
## ph.ecog 0.05947203 -0.19116847 0.23805821 0.30865378 -0.005363288
## ph.karno -0.02252266 0.09487913 -0.16127595 -0.32261297 -0.019623924
## pat.karno 0.04147893 0.17505701 -0.18542442 -0.23989736 0.071014942
## meal.cal 0.09869124 0.07467151 0.02483564 -0.23958240 -0.171044801
## wt.loss -0.17485406 0.03342528 0.04868879 0.04286056 -0.169892775
## ph.ecog ph.karno pat.karno meal.cal wt.loss
## inst 0.059472025 -0.02252266 0.04147893 0.09869124 -0.17485406
## time -0.191168469 0.09487913 0.17505701 0.07467151 0.03342528
## status 0.238058213 -0.16127595 -0.18542442 0.02483564 0.04868879
## age 0.308653782 -0.32261297 -0.23989736 -0.23958240 0.04286056
## sex -0.005363288 -0.01962392 0.07101494 -0.17104480 -0.16989278
## ph.ecog 1.000000000 -0.82269739 -0.54719617 -0.10563039 0.17125740
## ph.karno -0.822697393 1.00000000 0.53502749 0.05385409 -0.12524032
## pat.karno -0.547196168 0.53502749 1.00000000 0.17465190 -0.18213953
## meal.cal -0.105630385 0.05385409 0.17465190 1.00000000 -0.11134425
## wt.loss 0.171257402 -0.12524032 -0.18213953 -0.11134425 1.00000000#Escogemos las variables ph.ecog y ph.karno y hacemos el modelo de regresión lineal
regmod2<-lm(ph.ecog~ph.karno, data=data_lung)
summary(regmod2)##
## Call:
## lm(formula = ph.ecog ~ ph.karno, data = data_lung)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4658 -0.1127 -0.0539 0.4167 1.0049
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.818865 0.210176 22.93 <2e-16 ***
## ph.karno -0.047062 0.002532 -18.59 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4168 on 165 degrees of freedom
## Multiple R-squared: 0.6768, Adjusted R-squared: 0.6749
## F-statistic: 345.6 on 1 and 165 DF, p-value: < 2.2e-16lm(formula = ph.ecog ~ ph.karno, data = data_lung)##
## Call:
## lm(formula = ph.ecog ~ ph.karno, data = data_lung)
##
## Coefficients:
## (Intercept) ph.karno
## 4.81886 -0.04706plot(data_lung$ph.karno,data_lung$ph.ecog, xlab="ph.ecog", ylab="ph.karno")
abline(regmod2)
#Ahora haremos el modelo de regresión lineal múltiple
mrm_lung<-lm(data_lung$ph.ecog~data_lung$inst+data_lung$time+data_lung$status+data_lung$age+data_lung$sex+data_lung$ph.karno+data_lung$pat.karno+data_lung$meal.cal+data_lung$wt.loss, data=data_lung)
summary(mrm_lung)##
## Call:
## lm(formula = data_lung$ph.ecog ~ data_lung$inst + data_lung$time +
## data_lung$status + data_lung$age + data_lung$sex + data_lung$ph.karno +
## data_lung$pat.karno + data_lung$meal.cal + data_lung$wt.loss,
## data = data_lung)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.07299 -0.21742 -0.01057 0.29690 0.93285
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.434e+00 4.549e-01 9.746 <2e-16 ***
## data_lung$inst 6.777e-03 3.918e-03 1.730 0.0856 .
## data_lung$time -3.171e-04 1.528e-04 -2.075 0.0396 *
## data_lung$status 1.555e-01 7.299e-02 2.130 0.0347 *
## data_lung$age 1.133e-03 3.753e-03 0.302 0.7631
## data_lung$sex 3.141e-02 6.902e-02 0.455 0.6497
## data_lung$ph.karno -4.135e-02 3.001e-03 -13.779 <2e-16 ***
## data_lung$pat.karno -5.260e-03 2.525e-03 -2.083 0.0389 *
## data_lung$meal.cal -6.334e-05 8.166e-05 -0.776 0.4391
## data_lung$wt.loss 3.907e-03 2.443e-03 1.599 0.1119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3983 on 157 degrees of freedom
## Multiple R-squared: 0.7192, Adjusted R-squared: 0.7031
## F-statistic: 44.68 on 9 and 157 DF, p-value: < 2.2e-16#Seleccionamos los mejores predictores posibles
step(object=mrm_lung, direction="both",trace=1)## Start: AIC=-297.76
## data_lung$ph.ecog ~ data_lung$inst + data_lung$time + data_lung$status +
## data_lung$age + data_lung$sex + data_lung$ph.karno + data_lung$pat.karno +
## data_lung$meal.cal + data_lung$wt.loss
##
## Df Sum of Sq RSS AIC
## - data_lung$age 1 0.0145 24.924 -299.66
## - data_lung$sex 1 0.0329 24.942 -299.54
## - data_lung$meal.cal 1 0.0955 25.005 -299.12
## <none> 24.910 -297.76
## - data_lung$wt.loss 1 0.4056 25.315 -297.06
## - data_lung$inst 1 0.4747 25.384 -296.61
## - data_lung$time 1 0.6834 25.593 -295.24
## - data_lung$pat.karno 1 0.6883 25.598 -295.21
## - data_lung$status 1 0.7201 25.630 -295.00
## - data_lung$ph.karno 1 30.1251 55.035 -167.38
##
## Step: AIC=-299.66
## data_lung$ph.ecog ~ data_lung$inst + data_lung$time + data_lung$status +
## data_lung$sex + data_lung$ph.karno + data_lung$pat.karno +
## data_lung$meal.cal + data_lung$wt.loss
##
## Df Sum of Sq RSS AIC
## - data_lung$sex 1 0.026 24.951 -301.48
## - data_lung$meal.cal 1 0.126 25.050 -300.82
## <none> 24.924 -299.66
## - data_lung$wt.loss 1 0.399 25.323 -299.01
## - data_lung$inst 1 0.495 25.419 -298.38
## + data_lung$age 1 0.014 24.910 -297.76
## - data_lung$time 1 0.682 25.606 -297.15
## - data_lung$pat.karno 1 0.693 25.617 -297.08
## - data_lung$status 1 0.747 25.671 -296.73
## - data_lung$ph.karno 1 32.688 57.612 -161.73
##
## Step: AIC=-301.48
## data_lung$ph.ecog ~ data_lung$inst + data_lung$time + data_lung$status +
## data_lung$ph.karno + data_lung$pat.karno + data_lung$meal.cal +
## data_lung$wt.loss
##
## Df Sum of Sq RSS AIC
## - data_lung$meal.cal 1 0.158 25.109 -302.43
## <none> 24.951 -301.48
## - data_lung$wt.loss 1 0.375 25.326 -300.99
## - data_lung$inst 1 0.506 25.456 -300.13
## + data_lung$sex 1 0.026 24.924 -299.66
## + data_lung$age 1 0.008 24.942 -299.54
## - data_lung$time 1 0.661 25.612 -299.11
## - data_lung$pat.karno 1 0.677 25.628 -299.01
## - data_lung$status 1 0.720 25.671 -298.73
## - data_lung$ph.karno 1 33.216 58.166 -162.13
##
## Step: AIC=-302.43
## data_lung$ph.ecog ~ data_lung$inst + data_lung$time + data_lung$status +
## data_lung$ph.karno + data_lung$pat.karno + data_lung$wt.loss
##
## Df Sum of Sq RSS AIC
## <none> 25.109 -302.43
## - data_lung$wt.loss 1 0.414 25.523 -301.70
## + data_lung$meal.cal 1 0.158 24.951 -301.48
## - data_lung$inst 1 0.462 25.571 -301.38
## + data_lung$sex 1 0.059 25.050 -300.82
## + data_lung$age 1 0.033 25.076 -300.65
## - data_lung$status 1 0.674 25.783 -300.00
## - data_lung$time 1 0.703 25.812 -299.82
## - data_lung$pat.karno 1 0.802 25.911 -299.17
## - data_lung$ph.karno 1 33.085 58.194 -164.05##
## Call:
## lm(formula = data_lung$ph.ecog ~ data_lung$inst + data_lung$time +
## data_lung$status + data_lung$ph.karno + data_lung$pat.karno +
## data_lung$wt.loss, data = data_lung)
##
## Coefficients:
## (Intercept) data_lung$inst data_lung$time
## 4.5535875 0.0066252 -0.0003192
## data_lung$status data_lung$ph.karno data_lung$pat.karno
## 0.1466834 -0.0415905 -0.0055947
## data_lung$wt.loss
## 0.0038653#Ejecutamos de nuevo el modelo con los predictores más significativos
mrm_lung_predictores<-(lm(formula = data_lung$ph.ecog ~ data_lung$inst + data_lung$time +
data_lung$status + data_lung$ph.karno + data_lung$pat.karno +
data_lung$wt.loss, data = data_lung))
summary(mrm_lung_predictores)##
## Call:
## lm(formula = data_lung$ph.ecog ~ data_lung$inst + data_lung$time +
## data_lung$status + data_lung$ph.karno + data_lung$pat.karno +
## data_lung$wt.loss, data = data_lung)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.16209 -0.22765 -0.01216 0.30052 0.92432
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.5535875 0.2816769 16.166 <2e-16 ***
## data_lung$inst 0.0066252 0.0038613 1.716 0.0881 .
## data_lung$time -0.0003192 0.0001509 -2.116 0.0359 *
## data_lung$status 0.1466834 0.0707557 2.073 0.0398 *
## data_lung$ph.karno -0.0415905 0.0028644 -14.520 <2e-16 ***
## data_lung$pat.karno -0.0055947 0.0024744 -2.261 0.0251 *
## data_lung$wt.loss 0.0038653 0.0023803 1.624 0.1064
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3961 on 160 degrees of freedom
## Multiple R-squared: 0.7169, Adjusted R-squared: 0.7063
## F-statistic: 67.54 on 6 and 160 DF, p-value: < 2.2e-16#Realizamos un ANOVA
resultado_anova<-aov(sex ~ time, data = data_lung)
#Resumen del ANOVA
summary(resultado_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## time 1 0.51 0.5143 2.178 0.142
## Residuals 165 38.96 0.2361#Haremos un aAgrupamiento jerárquico divisivo. Primero calculamos la matriz de distancias
dist_matrix<-dist(data_lung)
#Cargamos el paquete "cluster" u usamos la función diana() para usar el método de agrupamiento divisivo
library(cluster)
hc_divisivo<-diana(dist_matrix) # Método de agrupamiento divisivo
# Representación del dendrograma
plot(hc_divisivo, main = "Dendrograma de Agrupamiento Jerárquico Divisivo", xlab = "data_lung")
#Elegimos un nivel de corte del dendograma, por ejemplo, dividir en 5 clústeres, como antes
nivel_corte<-4
cluster_asignado<-cutree(hc_divisivo, k = nivel_corte)
#Mostramos los resultados
resultado<-data.frame(Cluster = cluster_asignado)
print(resultado)## Cluster
## 1 1
## 2 1
## 3 2
## 4 3
## 5 3
## 6 3
## 7 3
## 8 1
## 9 4
## 10 2
## 11 2
## 12 3
## 13 1
## 14 1
## 15 1
## 16 2
## 17 2
## 18 1
## 19 1
## 20 3
## 21 1
## 22 3
## 23 1
## 24 1
## 25 2
## 26 2
## 27 1
## 28 2
## 29 1
## 30 1
## 31 3
## 32 2
## 33 3
## 34 3
## 35 2
## 36 2
## 37 2
## 38 1
## 39 1
## 40 2
## 41 3
## 42 2
## 43 1
## 44 3
## 45 1
## 46 1
## 47 1
## 48 3
## 49 2
## 50 1
## 51 3
## 52 2
## 53 3
## 54 1
## 55 3
## 56 1
## 57 2
## 58 1
## 59 3
## 60 2
## 61 1
## 62 1
## 63 3
## 64 2
## 65 1
## 66 2
## 67 1
## 68 2
## 69 1
## 70 2
## 71 1
## 72 3
## 73 3
## 74 2
## 75 2
## 76 1
## 77 1
## 78 2
## 79 1
## 80 1
## 81 3
## 82 1
## 83 1
## 84 3
## 85 3
## 86 1
## 87 3
## 88 3
## 89 3
## 90 2
## 91 3
## 92 2
## 93 1
## 94 1
## 95 1
## 96 2
## 97 1
## 98 1
## 99 3
## 100 3
## 101 1
## 102 1
## 103 1
## 104 3
## 105 1
## 106 1
## 107 3
## 108 3
## 109 1
## 110 3
## 111 1
## 112 1
## 113 3
## 114 3
## 115 1
## 116 3
## 117 4
## 118 4
## 119 1
## 120 3
## 121 1
## 122 1
## 123 3
## 124 1
## 125 1
## 126 1
## 127 3
## 128 3
## 129 1
## 130 1
## 131 1
## 132 3
## 133 3
## 134 1
## 135 1
## 136 1
## 137 1
## 138 3
## 139 1
## 140 1
## 141 1
## 142 3
## 143 3
## 144 1
## 145 1
## 146 3
## 147 1
## 148 1
## 149 3
## 150 4
## 151 1
## 152 3
## 153 3
## 154 1
## 155 3
## 156 3
## 157 1
## 158 1
## 159 3
## 160 3
## 161 3
## 162 1
## 163 1
## 164 4
## 165 1
## 166 1
## 167 1#Finalmente podemos añadir una representación gráfica de los diferentes clústers
pltree(hc_divisivo,hang=-1,cex=0.6)
rect.hclust(hc_divisivo,k=10,border=2:10)
#Caso práctico 3. ANOVA A partir del conjunto de datos birthwt del paquete MASS, vamos a intentar obtener conclusiones acerca de la relación existente entre la raza de la madre, teniendo en cuenta si esta es o no fumadora, y el bajo peso del recién nacido. Este estudio podría contextualizarse en una prueba ANOVA de una vía para datos independientes.
#Cargamos los datos
library(MASS)
data("birthwt")
head(birthwt)## low age lwt race smoke ptl ht ui ftv bwt
## 85 0 19 182 2 0 0 0 1 0 2523
## 86 0 33 155 3 0 0 0 0 3 2551
## 87 0 20 105 1 1 0 0 0 1 2557
## 88 0 21 108 1 1 0 0 1 2 2594
## 89 0 18 107 1 1 0 0 1 0 2600
## 91 0 21 124 3 0 0 0 0 0 2622#Como nos interesa estudiar el parámetro raza, cambiaremos los valores 1-2-3 por el nombre de la raza
birthwt$race<-factor(birthwt$race,levels=c(1,2,3),labels=c("white","black","other"))
#Ordenamos los datos que queremos usar para la ANOVA
raza<-c(birthwt$race)
peso<-c(birthwt$bwt)
df_birth_anova<-data.frame(raza,peso)
#Calculamos media y varianza
aggregate(peso~raza,data=df_birth_anova,FUN=mean)## raza peso
## 1 white 3102.719
## 2 black 2719.692
## 3 other 2805.284aggregate(peso~raza,data=df_birth_anova,FUN=sd)## raza peso
## 1 white 727.8861
## 2 black 638.6839
## 3 other 722.1944#Test para determinar la normalidad
require(nortest)## Loading required package: nortestby(data=df_birth_anova,INDICES=df_birth_anova$raza,FUN=function(x){lillie.test(x$peso)})## df_birth_anova$raza: white
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$peso
## D = 0.090389, p-value = 0.05113
##
## ------------------------------------------------------------
## df_birth_anova$raza: black
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$peso
## D = 0.1231, p-value = 0.3946
##
## ------------------------------------------------------------
## df_birth_anova$raza: other
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$peso
## D = 0.11598, p-value = 0.02578#Se podría asumir normalidad para white y black, pero no para others porque p-value=0.02578<0.05
#Test para determinar la varianza
fligner.test(df_birth_anova$peso~df_birth_anova$raza,df_birth_anova)##
## Fligner-Killeen test of homogeneity of variances
##
## data: df_birth_anova$peso by df_birth_anova$raza
## Fligner-Killeen:med chi-squared = 1.0086, df = 2, p-value = 0.6039#Hay homocedasticidad porque p-value=.6039>0.05
#Realizamos un ANOVA
resultados_birthwt_anova<-aov(df_birth_anova$peso~df_birth_anova$raza, data = df_birth_anova)
#Resumen del ANOVA
summary(resultados_birthwt_anova)## Df Sum Sq Mean Sq F value Pr(>F)
## df_birth_anova$raza 2 5015725 2507863 4.913 0.00834 **
## Residuals 186 94953931 510505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Gráficos de la ANOVA
plot(resultados_birthwt_anova)
#Realizamos un test de significancia
TukeyHSD(resultados_birthwt_anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = df_birth_anova$peso ~ df_birth_anova$raza, data = df_birth_anova)
##
## $`df_birth_anova$raza`
## diff lwr upr p adj
## black-white -383.02644 -756.2363 -9.816581 0.0428037
## other-white -297.43517 -566.1652 -28.705095 0.0260124
## other-black 85.59127 -304.4521 475.634630 0.8624372#Observamos que hay diferencias significativas entre black y white, y entre other y white; pero no entre other y black.