LAB5 - Introducció al Machine Learning amb R
Carles Castellà Barnaus
2026-05-14
EXERCICIS
EXERCICI 1
A partir del conjunt de dades swiss del paquet datasets, calculeu la matriu de correlació i, a partir d’aquesta, feu una anàlisi de regressió múltiple. Raoneu els resultats.
library(datasets)
data("swiss")
# Fem la matriu de correlació, per veure si hi ha alguna correlació alta entre variants
cor(swiss[,c("Fertility", "Agriculture", "Examination", "Education", "Catholic", "Infant.Mortality")], use ="complete", method = "pearson")## Fertility Agriculture Examination Education Catholic
## Fertility 1.0000000 0.35307918 -0.6458827 -0.66378886 0.4636847
## Agriculture 0.3530792 1.00000000 -0.6865422 -0.63952252 0.4010951
## Examination -0.6458827 -0.68654221 1.0000000 0.69841530 -0.5727418
## Education -0.6637889 -0.63952252 0.6984153 1.00000000 -0.1538589
## Catholic 0.4636847 0.40109505 -0.5727418 -0.15385892 1.0000000
## Infant.Mortality 0.4165560 -0.06085861 -0.1140216 -0.09932185 0.1754959
## Infant.Mortality
## Fertility 0.41655603
## Agriculture -0.06085861
## Examination -0.11402160
## Education -0.09932185
## Catholic 0.17549591
## Infant.Mortality 1.00000000# Al no detectar cap correlació alta, podem fer un anàlisi de regressió per acabar-ho de confirmar
model <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data = swiss)
summary(model)##
## Call:
## lm(formula = Fertility ~ Agriculture + Examination + Education +
## Catholic + Infant.Mortality, data = swiss)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2743 -5.2617 0.5032 4.1198 15.3213
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 66.91518 10.70604 6.250 1.91e-07 ***
## Agriculture -0.17211 0.07030 -2.448 0.01873 *
## Examination -0.25801 0.25388 -1.016 0.31546
## Education -0.87094 0.18303 -4.758 2.43e-05 ***
## Catholic 0.10412 0.03526 2.953 0.00519 **
## Infant.Mortality 1.07705 0.38172 2.822 0.00734 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared: 0.7067, Adjusted R-squared: 0.671
## F-statistic: 19.76 on 5 and 41 DF, p-value: 5.594e-10Com s’ha vist, primer en la matriu de correlació, cap variable té una correlació molt forta amb una altre, els valors més destacats son Education - Examination, amb un 0.69, i Agriculture - Examination, amb -0.68.
Els resultats del model de regressió lineal ens indiquen, primerament, que només un 70% de la variabilitat de la fertilitat s’explica pel model (R-squared), pel que hi ha variables incloses en el model que no son rellevants per l’estudi. Addicionalment, els valors de p-value indiquen que certes variables tenen significació en l’estudi, com education (la variable més influent, demostrat pel seu p-value), infant mortality i catholic (ambdues bastant significatives, però amb un efecte poc marcat (0.1 positiu)) i agriculture, que no és la més significativa. Per altra banda, la variable examination té un p-value molt superior a 0,05, pel que al ja saber el nivell d’educació (Education), la nota de l’examen militar no aporta informació extra, pel que es podria prescindir d’aquesta.
EXERCICI 2
A partir del conjunt de dades Pima.tr del paquet MASS que fa referència a dades de diabetis de pacients índies, feu un model de classificació SVM i feu una classificació de pacients per predir si estan en risc de patir diabetis. Per fer aquest exercici utilitzeu el paquet kernlab per entrenar el model SVM.
library(kernlab)
library(MASS)
data(Pima.tr)
set.seed(909)
# Dividim el dataset en dues parts
divisio <- sample(1:nrow(Pima.tr), 0.5*nrow(Pima.tr))
data_train <- Pima.tr[divisio, ]
data_test <- Pima.tr[- divisio, ]
# Definim el model a partir de les dades d'entrenament
model <- ksvm(type ~ ., data = data_train, kernel = "rbfdot")
# I la predicció a partir del model, sobre les dades de prova
predicts <- predict(model, newdata= data_test)
m_confusio <- table(Real = data_test$type, prediccio = predicts)
precisio <- sum(diag(m_confusio))/ sum(m_confusio)
cat("La precisió del model SVM és: ", precisio, "\n")## La precisió del model SVM és: 0.75m_confusio## prediccio
## Real No Yes
## No 55 17
## Yes 8 20EXERCICI 3
A partir del conjunt de dades Iris, es demana fer una anàlisi de components centrat, per exemple, en les mesures de longitud i amplada del sèpal i el pètal.
data(iris)
# Al tenir aquestes 4 categories numèriques, i la que no volem seleccionar que és tipus text, seleccionem amb is.numeric i estandaritzem al mateix temps
dades_iris <- scale(iris[, sapply(iris, is.numeric)])
head(dades_iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width
## [1,] -0.8976739 1.01560199 -1.335752 -1.311052
## [2,] -1.1392005 -0.13153881 -1.335752 -1.311052
## [3,] -1.3807271 0.32731751 -1.392399 -1.311052
## [4,] -1.5014904 0.09788935 -1.279104 -1.311052
## [5,] -1.0184372 1.24503015 -1.335752 -1.311052
## [6,] -0.5353840 1.93331463 -1.165809 -1.048667# Apliquem el PCA
pca_iris <- prcomp(dades_iris, center = TRUE, scale = TRUE)
summary(pca_iris)## Importance of components:
## PC1 PC2 PC3 PC4
## Standard deviation 1.7084 0.9560 0.38309 0.14393
## Proportion of Variance 0.7296 0.2285 0.03669 0.00518
## Cumulative Proportion 0.7296 0.9581 0.99482 1.00000plot(pca_iris)
Podem observar que PC1, és a dir, el component 1 recull més del 70% de la variabilitat total, mentre que PC2 recull un 22%. Entre aquests dos components ja s’explica el 95% de la informació de la flor, pel que les altres dues dimensions es podrien eliminar sense perdre massa informació (un 5%). Per tant, el dataset iris tot i ser quadridimensional, es pot simplificar a un pla bidimensional sense perdre la robustesa.
Per tant, podem concloure que les quatre variables estan molt correlacionades entre elles (de no ser així, al ser independents entre elles cada component recolliria un 25% de la variança).
EXERCICI 4
Feu un agrupament jeràrquic aglomeratiu amb el paquet de dades birthwt del paquet MASS.
library(cluster)
library(MASS)
data("birthwt")
set.seed(999)
# Estandaritzem les dades, i guardem només les dades numèriques
dades <- scale(birthwt[, sapply(birthwt, is.numeric)])
# Calculem la distància entre cada parell de pacients
distancies <- dist(dades, method = "euclidean")
# Fem l'agrupament jeràrquic mitjançant Ward
hc_model <- hclust(distancies, method = "ward.D2")
# I representem gràficament
plot(hc_model, main = "Dendrograma d'Agrupament Jeràrquic Aglomeratiu", xlab = "Pacients")
# També podem calcular-ho amb la funció agnes, per obtenir el coeficient;
hc_agnes <- agnes(dades, method = "ward")
coef_aglom <- hc_agnes$ac
cat("El coeficient d'aglomeració és: ", coef_aglom, "\n")## El coeficient d'aglomeració és: 0.9387985Obtenim que el coeficient d’agloimeració és molt proper a 1, pel que la fortalesa de l’agrupament és molt alta, és a dir, la distància entre els individus i els clústers està ben definida, validant que aquesta estratègia aglomerativa és adequada.
EXERCICI 5
Feu un agrupament jeràrquic divisiu amb el paquet de dades birthwt del paquet MASS.
library(cluster)
library(MASS)
data("birthwt")
set.seed(606)
# Estandaritzem les dades, i guardem només les dades numèriques
dades <- scale(birthwt[, sapply(birthwt, is.numeric)])
hc_divisiu <- diana(dades)
pltree(hc_divisiu,cex=0.5,hang=-1,main="Dendrograma de Diana")
coef_div <- hc_divisiu$dc
cat("El coeficient de divisió és: ", coef_div, "\n")## El coeficient de divisió és: 0.8481872De nou, hem obtingut un coeficient bastant proper a 1, tot i que en aquesta ocasió una mica inferior. Tot i així, és un bon valor, pel que se suposa una estructura divisiva forta, pel que ha estat fàcil per DIANA anar separant la mostra en grups.
EXERCICI 6
Feu un agrupament jeràrquic de tipus aglomeratiu per al conjunt de dades melanoma del paquet MASS.
library(cluster)
library(MASS)
data("Melanoma")
set.seed(988)
# Estandaritzem les dades, i guardem només les dades numèriques
dades <- scale(Melanoma[, c("time", "age", "thickness")])
# Calculem la distància entre cada parell de pacients
distancies <- dist(dades, method = "euclidean")
# Fem l'agrupament jeràrquic mitjançant Ward
hc_model <- hclust(distancies, method = "ward.D2")
# I representem gràficament
plot(hc_model, main = "Dendrograma d'Agrupament Jeràrquic Aglomeratiu", xlab = "Pacients")
# També podem calcular-ho amb la funció agnes, per obtenir el coeficient;
hc_agnes <- agnes(dades, method = "ward")
coef_aglom <- hc_agnes$ac
cat("El coeficient d'aglomeració és: ", coef_aglom, "\n")## El coeficient d'aglomeració és: 0.9771963Obtenim que el coeficient d’agloimeració és molt proper a 1, pel que la fortalesa de l’agrupament és molt alta, és a dir, la distància entre els individus i els clústers està ben definida, validant que aquesta estratègia aglomerativa és adequada.
EXERCICI 7
Feu un agrupament jeràrquic de tipus divisiu per al conjunt de dades melanoma del paquet MASS.
library(cluster)
library(MASS)
data("Melanoma")
set.seed(696)
# Estandaritzem les dades, i guardem només les dades numèriques
dades <- scale(Melanoma[, c("time", "age", "thickness")])
hc_divisiu <- diana(dades)
pltree(hc_divisiu,cex=0.5,hang=-1,main="Dendrograma de Diana")
coef_div <- hc_divisiu$dc
cat("El coeficient de divisió és: ", coef_div, "\n")## El coeficient de divisió és: 0.9308206EXERCICI 8
Feu un agrupament no jeràrquic utilitzant el conjunt de dades birthwt del paquet MASS.
library(MASS)
library(cluster)
data("birthwt")
dades <- MASS::birthwt
pairs(dades[, c("age", "lwt", "bwt")], col = dades$smoke)
corr <- cor(dades)
print(corr)## low age lwt race smoke ptl
## low 1.00000000 -0.11893928 -0.16962694 0.137792751 0.16140431 0.196087267
## age -0.11893928 1.00000000 0.18007315 -0.172817953 -0.04434618 0.071606386
## lwt -0.16962694 0.18007315 1.00000000 -0.165048544 -0.04417908 -0.140029003
## race 0.13779275 -0.17281795 -0.16504854 1.000000000 -0.33903074 0.007951293
## smoke 0.16140431 -0.04434618 -0.04417908 -0.339030745 1.00000000 0.187557063
## ptl 0.19608727 0.07160639 -0.14002900 0.007951293 0.18755706 1.000000000
## ht 0.15237025 -0.01583700 0.23636040 0.019929917 0.01340704 -0.015399579
## ui 0.16904283 -0.07515558 -0.15276317 0.053602088 0.06215900 0.227585340
## ftv -0.06296026 0.21539394 0.14052746 -0.098336254 -0.02801314 -0.044429660
## bwt -0.78480516 0.09031781 0.18573328 -0.194713487 -0.19044806 -0.154653390
## ht ui ftv bwt
## low 0.15237025 0.16904283 -0.06296026 -0.78480516
## age -0.01583700 -0.07515558 0.21539394 0.09031781
## lwt 0.23636040 -0.15276317 0.14052746 0.18573328
## race 0.01992992 0.05360209 -0.09833625 -0.19471349
## smoke 0.01340704 0.06215900 -0.02801314 -0.19044806
## ptl -0.01539958 0.22758534 -0.04442966 -0.15465339
## ht 1.00000000 -0.10858506 -0.07237255 -0.14598189
## ui -0.10858506 1.00000000 -0.05952341 -0.28392741
## ftv -0.07237255 -0.05952341 1.00000000 0.05831777
## bwt -0.14598189 -0.28392741 0.05831777 1.00000000# Escalem les dades per estandaritzar els valors i evitar que s'emmascarin resultats
scaled <- scale(dades)
set.seed(586)
grups_km <- kmeans(scaled, centers = 2)
grups_km$cluster## 85 86 87 88 89 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
## 2 2 2 1 1 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2
## 106 107 108 109 111 112 113 114 115 116 117 118 119 120 121 123 124 125 126 127
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
## 2 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 148 149 150 151 154 155 156 159 160 161 162 163 164 166 167 168 169 170 172 173
## 2 2 2 2 1 2 2 2 1 2 1 2 2 2 2 2 2 2 2 2
## 174 175 176 177 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 195
## 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2
## 196 197 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 217 218 219 220 221 222 223 224 225 226 4 10 11 13 15 16 17 18 19 20
## 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1
## 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 40 42 43 44
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 45 46 47 49 50 51 52 54 56 57 59 60 61 62 63 65 67 68 69 71
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 75 76 77 78 79 81 82 83 84
## 1 1 1 1 1 1 1 1 1plot(grups_km$cluster)
table(Fumadora = birthwt$smoke, Cluster = grups_km$cluster)## Cluster
## Fumadora 1 2
## 0 31 84
## 1 37 37L’anàlisi de clúster K-means (K = 2) realitzada sobre les dades estandarditzades de birthwt ens permet identificar dos perfils de pacients. En comparar aquests clústers amb la variable smoke, observem si l’estil de vida (tabac) és un factor determinant en la configuració bionatal o si, per contra, la variabilitat del conjunt de dades està dirigida per altres variables biomètriques amb més pes estadístic. Veiem en executar l’ordre corr, que hi ha una correlació molt baixa entre bwt i smoke, mentre que aquesta és molt alta entre low (pes baix) i bwt (pes en grams), el qual té sentit perque una deriva de l’altre.
Trobem també (resultats de la taula) que l’algoritme no ha sigut capaç de reconstruir el grup de fumadores (resultats no disposats en diagonal), pel que el perfil de les fumadores no es diferencia de les que no ho son, pel que es pot concloure que el fet de fumar no és el factor principal que defineix la similitud o les diferències entre pacients d’aquest dataset.
EXERCICI 9
Feu un agrupament no jeràrquic utilitzant la tècnica k-means() per al conjunt de dades melanoma del paquet MASS. Valoreu els resultats per a diferents nombres de clústers i definint una determinada variable objectiu. Podeu transformar alguna variable sencera de manera categòrica, per exemple, la referent a l’estatus.
library(MASS)
library(cluster)
data("Melanoma")
# Mirem les diferents variables que hi ha
head(Melanoma)## time status sex age year thickness ulcer
## 1 10 3 1 76 1972 6.76 1
## 2 30 3 1 56 1968 0.65 0
## 3 35 2 1 41 1977 1.34 0
## 4 99 3 0 71 1968 2.90 0
## 5 185 1 1 52 1965 12.08 1
## 6 204 1 1 28 1971 4.84 1dades <- Melanoma
# Fem un pairs excloent la variable status, doncs tot i ser numèrica, es categòrica. # Llavors dividim els resultats segons l'estat del pacient, així podem veure si els punts d'un mateix estat es concentren en algun lloc
pairs(dades[-2],col=Melanoma$status)
# Comprovem si hi ha alguna correlació entre les diferents variables
corr <- cor(Melanoma[-2])
corr## time sex age year thickness
## time 1.0000000 -0.146499215 -0.30151794 -0.485504359 -0.2354087
## sex -0.1464992 1.000000000 0.06833741 -0.002645159 0.1854126
## age -0.3015179 0.068337413 1.00000000 0.188229089 0.2124798
## year -0.4855044 -0.002645159 0.18822909 1.000000000 -0.1333454
## thickness -0.2354087 0.185412563 0.21247979 -0.133345424 1.0000000
## ulcer -0.2647575 0.167979154 0.12606294 -0.033125618 0.4244593
## ulcer
## time -0.26475748
## sex 0.16797915
## age 0.12606294
## year -0.03312562
## thickness 0.42445931
## ulcer 1.00000000# Com es pot observar, hi ha una correlació negativa entre time i thickness/ulcer. El sentit clínic d'això és que a major gruix del tumor (o si aquest s'ulcera, que té una relació positiva de 0.42 amb thickness), menor és el temps de supervivència del pacient.
# A partir d'això, trobem els k-means
set.seed(654)
grups <- kmeans(Melanoma[-2], centers = 3)
grups$cluster## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 3 3 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 1
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 201 202 203 204 205
## 1 1 1 1 1plot(grups$cluster)
table(Estatus_real = dades$status, Cluster_trobat = grups$cluster)## Cluster_trobat
## Estatus_real 1 2 3
## 1 2 35 20
## 2 50 1 83
## 3 3 7 4L’algoritme és capaç de fer una predicció o segmentació no supervisada. L’alta correspondència que mostra la taula final entre els clústers i l’estatus real del pacient valida que les variables físiques i de diagnòstic (thickness, time, ulcer) són suficients per descriure de manera inherent el patró de gravetat del melanoma, fins i tot sense conèixer el desenllaç de l’estudi.
EXERCICI 10
A partir del conjunt de dades Iris, apliqueu la funció ANOVA i raoneu els resultats obtinguts.
data("iris")
head(iris)## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 5.1 3.5 1.4 0.2 setosa
## 2 4.9 3.0 1.4 0.2 setosa
## 3 4.7 3.2 1.3 0.2 setosa
## 4 4.6 3.1 1.5 0.2 setosa
## 5 5.0 3.6 1.4 0.2 setosa
## 6 5.4 3.9 1.7 0.4 setosaanova <- aov(Sepal.Width ~ Species, data = iris)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## Species 2 11.35 5.672 49.16 <2e-16 ***
## Residuals 147 16.96 0.115
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey <- TukeyHSD(anova)
tukey## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Sepal.Width ~ Species, data = iris)
##
## $Species
## diff lwr upr p adj
## versicolor-setosa -0.658 -0.81885528 -0.4971447 0.0000000
## virginica-setosa -0.454 -0.61485528 -0.2931447 0.0000000
## virginica-versicolor 0.204 0.04314472 0.3648553 0.0087802El p-valor obtingut és inferior al nivell de significació 0.05, pel que es pot rebutjar la hipòtesi nul·la d’igualtat de mitjanes, és a dir, hi ha diferències significatives en l’amplada del sepal entre almenys dues de les espècies, i no son produïdes per l’atzar. Posteriorment, he fet un test de Tukey que ha confirmat que les diferències existeixen i son presents en tots els parells d’espècies comparats
CASOS PRÀCTICS
CAS PRÀCTIC 1
A partir del conjunt de dades Ovarian del paquet survival, feu un estudi de regressió lineal sobre un parell de variables que, a priori, penseu que poden presentar correlació, i un estudi de regressió múltiple seleccionant els millors predictors.
library(survival)
data(ovarian)## Warning in data(ovarian): data set 'ovarian' not found# Fem un head per mostrar quines variables hi ha
head(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 2# Fem un gràfic pairs per veure si intuim alguna relació:
pairs(~ovarian$futime + ovarian$fustat + ovarian$age + ovarian$resid.ds + ovarian$rx + ovarian$ecog.ps)
# L'únic que no està segmentat i presenta una mica de relació amb pendent és age i futime, pel que n'analitzem la correlació
corr <- cor(ovarian[,c("age", "futime")], use = "complete")
print(corr)## age futime
## age 1.0000000 -0.6483612
## futime -0.6483612 1.0000000# Veiem que hi ha una correlació negativa de -0.648, que tot i no ser molt propera a -1, és prou considerable.
# Per tant, farem el model a partir d'aquestes dues variables
model <- lm(futime~age, data = ovarian)
summary(model)##
## Call:
## lm(formula = futime ~ age, data = ovarian)
##
## Residuals:
## Min 1Q Median 3Q Max
## -554.63 -160.13 -49.27 67.21 702.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1824.240 298.079 6.120 2.54e-06 ***
## age -21.805 5.227 -4.172 0.000341 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 263.9 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.0003408# Obtenim un r-squared no massa alt, de 0.42, pel que el nostre model no està massa ajustat a les dades. Per altra banda, el p-value és molt inferior al valor significatiu de 0.05, pel que ambdues variables tenen una relació.
# Ara farem el model de regressió múltiple, on considerem totes les variables
# Primer estimem el coeficient de correlació de pearson per totes les 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# L'única que té un coeficient de correlació similar és fustat
# Construim el model múltiple
model_m <- lm(futime ~ fustat + age + resid.ds + rx + ecog.ps, data = ovarian)
summary(model_m)##
## Call:
## lm(formula = futime ~ fustat + age + resid.ds + rx + 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 ***
## fustat -309.207 105.682 -2.926 0.008357 **
## age -14.373 5.083 -2.828 0.010401 *
## resid.ds -48.136 95.331 -0.505 0.619126
## rx 125.638 87.254 1.440 0.165364
## 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# Ara veiem que l'R-squared ha pujat fins a 0.67, un augment significatiu que indica que el model ara està molt millor ajustat a les dades. El p-value, en aquest cas, segueix per sota del nivel 0.05.
# Seleccionem els millors predictors mitjançant step():
step(object = model_m, direction = "both", trace = 1)## Start: AIC=285
## futime ~ fustat + age + resid.ds + rx + ecog.ps
##
## Df Sum of Sq RSS AIC
## - resid.ds 1 12034 956070 283.32
## - ecog.ps 1 69130 1013166 284.83
## <none> 944036 285.00
## - rx 1 97865 1041901 285.56
## - age 1 377400 1321436 291.74
## - fustat 1 404069 1348105 292.26
##
## Step: AIC=283.32
## futime ~ fustat + age + rx + ecog.ps
##
## Df Sum of Sq RSS AIC
## <none> 956070 283.32
## - ecog.ps 1 88920 1044991 283.64
## - rx 1 101328 1057398 283.94
## + resid.ds 1 12034 944036 285.00
## - age 1 414281 1370351 290.69
## - fustat 1 467936 1424007 291.68##
## Call:
## lm(formula = futime ~ fustat + age + rx + ecog.ps, data = ovarian)
##
## Coefficients:
## (Intercept) fustat age rx ecog.ps
## 1213.26 -322.41 -14.82 127.70 120.52# Ara cridem la funció que ens ha especificat l'step()
millor_model <- lm(formula = futime ~ fustat + age + rx + ecog.ps, data = ovarian)
summary(millor_model)##
## Call:
## lm(formula = futime ~ fustat + age + rx + 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 ***
## fustat -322.407 100.565 -3.206 0.004244 **
## age -14.824 4.914 -3.017 0.006568 **
## rx 127.701 85.599 1.492 0.150606
## 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-05En primer lloc, l’estudi de regressió lineal simple revela una forta i significativa correlació lineal negativa entre l’edat i el temps de supervivència (r = -0.648, i p < 0.001).
En segon lloc, l’aplicació de l’algorisme de selecció per passes (step())ha optimitzat el model multivariant eliminant la variable de malaltia residual (resid.ds) que considera redundant. El model múltiple guanyador reté fustat, age, tx i ecog.ps assolint una capacitat explicativa amb R-quadrat ajustat del 60,54% i p < 0.001.
CAS PRÀCTIC 2
A partir del conjunt de dades lung del paquet survival, es demana fer un breu estudi aplicant els conceptes estudiats en aquest LAB. Per fer-ho, es demana dur a terme un estudi de regressió lineal, regressió múltiple, test ANOVA i test clustering sobre aquest conjunt de dades. Les variables triades en cadascun dels casos possibles a tractar es donen a triar a l’estudiant.
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# Veiem que hi ha valors buits (NA), pel que hem de netejar-ho tot per evitar errors
dades <- na.omit(lung)
# Busquem mitjançant pairs si hi ha alguna relació entre les variables
pairs(~inst + time + status + age + sex + ph.ecog + ph.karno + pat.karno + meal.cal + wt.loss, data=dades)
# Per representar les relacions numèricament, estimem la correlació
corr <- cor(x = dades, method = "pearson")
print(corr)## 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# Les variables que tenen una correlació més clara son ph.ecog i ph.karno, amb una correlació negativa -0.822. Pat.karno també mostra certa correlació amb ph.egoc (-0.547) i amb ph.karno (0.535).
# Farem el model amb els dos ph.
model <- lm(formula = ph.ecog ~ ph.karno, data = dades)
summary(model)##
## Call:
## lm(formula = ph.ecog ~ ph.karno, data = dades)
##
## 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-16# El model presenta un R-squared de 0.67, pel que té un bon ajust sobre les dades. Farem el model múltiple per veure si altres variables influeixen en l'ajust
model_m <- lm(ph.ecog ~ inst + time + status + age + sex + ph.karno + pat.karno + meal.cal + wt.loss, data= dades)
summary(model_m)##
## Call:
## lm(formula = ph.ecog ~ inst + time + status + age + sex + ph.karno +
## pat.karno + meal.cal + wt.loss, data = dades)
##
## 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 ***
## inst 6.777e-03 3.918e-03 1.730 0.0856 .
## time -3.171e-04 1.528e-04 -2.075 0.0396 *
## status 1.555e-01 7.299e-02 2.130 0.0347 *
## age 1.133e-03 3.753e-03 0.302 0.7631
## sex 3.141e-02 6.902e-02 0.455 0.6497
## ph.karno -4.135e-02 3.001e-03 -13.779 <2e-16 ***
## pat.karno -5.260e-03 2.525e-03 -2.083 0.0389 *
## meal.cal -6.334e-05 8.166e-05 -0.776 0.4391
## 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# Trobem els millors predictors
step(object = model_m, direction = "both", trace = 1)## Start: AIC=-297.76
## ph.ecog ~ inst + time + status + age + sex + ph.karno + pat.karno +
## meal.cal + wt.loss
##
## Df Sum of Sq RSS AIC
## - age 1 0.0145 24.924 -299.66
## - sex 1 0.0329 24.942 -299.54
## - meal.cal 1 0.0955 25.005 -299.12
## <none> 24.910 -297.76
## - wt.loss 1 0.4056 25.315 -297.06
## - inst 1 0.4747 25.384 -296.61
## - time 1 0.6834 25.593 -295.24
## - pat.karno 1 0.6883 25.598 -295.21
## - status 1 0.7201 25.630 -295.00
## - ph.karno 1 30.1251 55.035 -167.38
##
## Step: AIC=-299.66
## ph.ecog ~ inst + time + status + sex + ph.karno + pat.karno +
## meal.cal + wt.loss
##
## Df Sum of Sq RSS AIC
## - sex 1 0.026 24.951 -301.48
## - meal.cal 1 0.126 25.050 -300.82
## <none> 24.924 -299.66
## - wt.loss 1 0.399 25.323 -299.01
## - inst 1 0.495 25.419 -298.38
## + age 1 0.014 24.910 -297.76
## - time 1 0.682 25.606 -297.15
## - pat.karno 1 0.693 25.617 -297.08
## - status 1 0.747 25.671 -296.73
## - ph.karno 1 32.688 57.612 -161.73
##
## Step: AIC=-301.48
## ph.ecog ~ inst + time + status + ph.karno + pat.karno + meal.cal +
## wt.loss
##
## Df Sum of Sq RSS AIC
## - meal.cal 1 0.158 25.109 -302.43
## <none> 24.951 -301.48
## - wt.loss 1 0.375 25.326 -300.99
## - inst 1 0.506 25.456 -300.13
## + sex 1 0.026 24.924 -299.66
## + age 1 0.008 24.942 -299.54
## - time 1 0.661 25.612 -299.11
## - pat.karno 1 0.677 25.628 -299.01
## - status 1 0.720 25.671 -298.73
## - ph.karno 1 33.216 58.166 -162.13
##
## Step: AIC=-302.43
## ph.ecog ~ inst + time + status + ph.karno + pat.karno + wt.loss
##
## Df Sum of Sq RSS AIC
## <none> 25.109 -302.43
## - wt.loss 1 0.414 25.523 -301.70
## + meal.cal 1 0.158 24.951 -301.48
## - inst 1 0.462 25.571 -301.38
## + sex 1 0.059 25.050 -300.82
## + age 1 0.033 25.076 -300.65
## - status 1 0.674 25.783 -300.00
## - time 1 0.703 25.812 -299.82
## - pat.karno 1 0.802 25.911 -299.17
## - ph.karno 1 33.085 58.194 -164.05##
## Call:
## lm(formula = ph.ecog ~ inst + time + status + ph.karno + pat.karno +
## wt.loss, data = dades)
##
## Coefficients:
## (Intercept) inst time status ph.karno pat.karno
## 4.5535875 0.0066252 -0.0003192 0.1466834 -0.0415905 -0.0055947
## wt.loss
## 0.0038653# I apliquem la fòrmula que ens ha donat l'step()
millor_model <- lm(formula = ph.ecog ~ inst + time + status + ph.karno + pat.karno + wt.loss, data = dades)
summary(millor_model)##
## Call:
## lm(formula = ph.ecog ~ inst + time + status + ph.karno + pat.karno +
## wt.loss, data = dades)
##
## 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 ***
## inst 0.0066252 0.0038613 1.716 0.0881 .
## time -0.0003192 0.0001509 -2.116 0.0359 *
## status 0.1466834 0.0707557 2.073 0.0398 *
## ph.karno -0.0415905 0.0028644 -14.520 <2e-16 ***
## pat.karno -0.0055947 0.0024744 -2.261 0.0251 *
## 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# Ara obtenim un r-squared encara millor, de 0.7169, que explica millor l'ajust sobre les dades. Hem vist que els predictors més útils per aquest model son ph.ecog, inst, time, status, ph.karno, pat.karno i wt.loss.
# Ara farem un ANOVA sobre les dues variables amb major correlació:
# Convertim l'ECOG a factor (degut a que es categòrica)
dades$dades_ecog <- factor(dades$ph.ecog)
anova <- aov(ph.karno ~ dades_ecog, data = dades)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## dades_ecog 3 18553 6184 117.8 <2e-16 ***
## Residuals 163 8555 52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1TukeyHSD(anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ph.karno ~ dades_ecog, data = dades)
##
## $dades_ecog
## diff lwr upr p adj
## 1-0 -12.760704 -16.20890 -9.312506 0.0000000
## 2-0 -29.316909 -33.41940 -25.214421 0.0000000
## 3-0 -35.106383 -54.11061 -16.102153 0.0000215
## 2-1 -16.556205 -20.25379 -12.858624 0.0000000
## 3-1 -22.345679 -41.26663 -3.424726 0.0134375
## 3-2 -5.789474 -24.84053 13.261584 0.8594505# Finalment, fem l'agrupament jeràrquic (clustering)
library(MASS)
library(cluster)
#Calculem les distàncies
distancia <- dist(dades, method = "euclidean")
cluster <- hclust(distancia, method = "complete")
#I representem gràficament
plot(cluster, main = "Dendrograma d'Agrupament Jeràrquic Aglomeratiu", xlab = "Pacients")
# També podem calcular-ho amb la funció agnes, per obtenir el coeficient;
hc_agnes <- agnes(dades, method = "ward")
coef_aglom <- hc_agnes$ac
cat("El coeficient d'aglomeració és: ", coef_aglom, "\n")## El coeficient d'aglomeració és: 0.987616L’anàlisi conjunta de regressió i variància (ANOVA) confirma una
congruència clínica absoluta entre les dues principals escales de
valoració oncològica del dataset lung. El model de
regressió múltiple assoleix una gran capacitat explicativa (\(R^2 \text{ ajustat} = 0.7063, p <
0.001\)), on l’índex de Karnofsky del metge (ph.karno)
s’estableix com el predictor més robust de l’estat funcional.
Per la seva banda, l’ANOVA corrobora que la mitjana de l’índex de Karnofsky varia significativament en funció dels graus de l’escala ECOG (F = 117.8, p < 0.001). El test post-hoc de Tukey detalla un descens progressiu i altament significatiu de la puntuació de Karnofsky a mesura que s’agreuja l’ECOG
Tanmateix, s’identifica una excepció rellevant en el contrast entre els graus ECOG 2 i 3 (p = 0.8594), indicant que l’estat biològic de pèrdua d’autonomia severa es troba estancat en un llindar de Karnofsky similar per a ambdós subgrups de pacients.
CAS PRÀCTIC 3
A partir del conjunt de dades birthwt del paquet MASS, intentarem obtenir conclusions sobre la relació existent entre la raça de la mare, tenint en compte si aquesta és o no fumadora, i el baix pes del nounat. Aquest estudi podria contextualitzar-se en una prova ANOVA d’una via per a dades independents.
library(cluster)
library(MASS)
data("birthwt")
dades <- birthwt
#Factoritzem les dades de raça i fumadores (HE INTENTAT FER L'EXERCICI INCLOENT
# FUMADORES, PERÒ QUAN HE VOLGUT CONSULTAR SI ESTAVA BÉ, AL SOLUCIONARI NO ES TÉ EN # COMPTE, PEL QUE HO HE FET COM AL SOLUCIONARI)
dades$race <- factor(dades$race, levels = c(1, 2, 3),
labels = c("Blanca", "Negra", "Altra"))
race <- c(dades$race)
pes <- c(dades$bwt)
# Definim les dades com un dataframe
df_birth_anova <- data.frame(race, pes)
table(df_birth_anova$race)##
## Blanca Negra Altra
## 96 26 67# Calculem mitjana i desviació segons race
tapply(df_birth_anova$pes, df_birth_anova$race, mean)## Blanca Negra Altra
## 3102.719 2719.692 2805.284tapply(df_birth_anova$pes, df_birth_anova$race, sd)## Blanca Negra Altra
## 727.8861 638.6839 722.1944# Abans d'aplicar un ANOVA s'han de complir dues condicions: la mostra ha de ser normal, i hi ha d'haver homoscedasticitat
# NORMALITAT:
require(nortest)## Cargando paquete requerido: nortestby(data=df_birth_anova, INDICES = df_birth_anova$race, FUN=function(x){lillie.test(x$pes)})## df_birth_anova$race: Blanca
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$pes
## D = 0.090389, p-value = 0.05113
##
## ------------------------------------------------------------
## df_birth_anova$race: Negra
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$pes
## D = 0.1231, p-value = 0.3946
##
## ------------------------------------------------------------
## df_birth_anova$race: Altra
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: x$pes
## D = 0.11598, p-value = 0.02578# HOMOSCEDASTICITAT
fligner.test(df_birth_anova$pes~df_birth_anova$race, df_birth_anova)##
## Fligner-Killeen test of homogeneity of variances
##
## data: df_birth_anova$pes by df_birth_anova$race
## Fligner-Killeen:med chi-squared = 1.0086, df = 2, p-value = 0.6039# ANOVA
anova <- aov(df_birth_anova$pes~df_birth_anova$race, data = df_birth_anova)
summary(anova)## Df Sum Sq Mean Sq F value Pr(>F)
## df_birth_anova$race 2 5015725 2507863 4.913 0.00834 **
## Residuals 186 94953931 510505
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(anova)
TukeyHSD(anova)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = df_birth_anova$pes ~ df_birth_anova$race, data = df_birth_anova)
##
## $`df_birth_anova$race`
## diff lwr upr p adj
## Negra-Blanca -383.02644 -756.2363 -9.816581 0.0428037
## Altra-Blanca -297.43517 -566.1652 -28.705095 0.0260124
## Altra-Negra 85.59127 -304.4521 475.634630 0.8624372Finalment, els resultats de l’ANOVA determinen que els fills de mares de raça negra tenen un pes significativament menor en nèixer que els de raça blanca, al igual que els de raça ‘altra’. No existeien diferències significatives entre el pes dels nadons de mares de raça negra i d’altres, ambdos es comporten de forma homogènia respecte a la reducció de pes.
Tot i que metodològicament s’ha seguit el criteri del solucionari executant un ANOVA d’una via centrat exclusivament en l’efecte aïllat de la raça, una anàlisi bivariant prèvia suggereix que fumar (smoke) actua com un factor rellevant. Per a futures anàlisis i per evitar un possible biaix de variable omesa, seria recomanable plantejar un ANOVA de dues vies amb interacció per comparar si la caiguda del pes en els grups de risc està més lligada a la biologia o als hàbits de salut de la mare.