Un model estadístic és una representació matemàtica que relaciona una variable dependent amb una o més variables independents. Té aquesta estructura general:
\[ Y = f(X) + \varepsilon\ \] On:
⚙️ 2. Tipus bàsics de models
🔹 Models lineals
Regressió lineal simple: \(Y= \beta_0 + \beta_1X + \varepsilon\)
Regressió lineal múltiple: \(Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \dots + \varepsilon\)
🔹 Models amb variables categòriques
🔹 Models logístics
🧪 3. Objectius d’un model
Imagina que volem predir l’alçada d’una persona, però no sabem res més d’ella. L’únic que tenim són les alçades d’un grup de persones. 💡 En aquest cas, el millor model possible (en absència de més informació) és la mitjana.
set.seed(42)
alçades <- rnorm(100, mean = 170, sd = 10)
mitjana <- mean(alçades)
plot(alçades, rep(1, 100), pch = 16, col = "gray", yaxt = "n",
xlab = "Alçada (cm)", ylab = "", main = "La mitjana com a model")
abline(v = mitjana, col = "blue", lwd = 2)
Aquest model ens diu: “Si no sé res més, et diré que tothom fa 170 cm”. Però no encertarem mai exactament: hi ha variabilitat.
Cada persona s’allunya de la mitjana. Aquesta diferència s’anomena residu. El conjunt dels residus ens indica com d’allunyat està el món real del nostre model.
residus <- alçades - mitjana
hist(residus, col = "lightblue", main = "Distribució dels residus",
xlab = "Alçada real - Mitjana", breaks = 10)
La desviació estàndard (SD) és una mesura del promig dels errors. Com més gran sigui la SD, pitjor és el nostre model.
Ara imagina que sabem una altra cosa: l’edat de cada persona. Podem fer un model millor? Sí. Ara ja no ens conformem amb una mitjana constant. Podem ajustar una recta que s’adapti als canvis d’edat:
set.seed(123)
edat <- sample(20:80, 100, replace = TRUE)
alçada <- 150 + 0.4 * edat + rnorm(100, 0, 5)
dades <- data.frame(edat, alçada)
model <- lm(alçada ~ edat, data = dades)
plot(dades$edat, dades$alçada, pch = 16, col = "darkgray",
xlab = "Edat (anys)", ylab = "Alçada (cm)",
main = "Model de regressió lineal")
abline(model, col = "blue", lwd = 2)
Ara el nostre model diu: “Segons l’edat, t’estimo una alçada”. Ja no és una línia recta horitzontal (la mitjana), sinó una recta amb pendent. I com abans, els punts no encaixen exactament amb la línia: hi ha residus.
En aquest cas analitzant la variabilitat tindríem:
Variabilitat total = quant varien les alçades (variança total)
Explicada: part que explica el model (canvis per edat)
No explicada: els residus
La mitjana és un model simple, però útil. La regressió afegeix predictors per fer millors estimacions.
Sempre hi ha error: la clau és quant n’expliquem i quant en queda. Un model millor té residus més petits i explica més del total.
Aprendrem a ajustar un model de regressió lineal simple, interpretar els coeficients, avaluar la seva significació i validar el model gràficament.
En la regressió lineal simple modelitzem la relació entre dues variables quantitatives:
\[ Y = \beta\_0 + \beta\_1 X + \varepsilon \]
On:
Analitzem la relació entre edat i pressió arterial sistòlica:
# Creem une dades d'entrenament
set.seed(123)
n <- 100
edat <- sample(30:80, n, replace = TRUE)
pressio <- 90 + 0.7 * edat + rnorm(n, 0, 10)
dades <- data.frame(edat, pressio)
# construim el model
model <- lm(pressio ~ edat, data = dades)
summary(model)
Call:
lm(formula = pressio ~ edat, data = dades)
Residuals:
Min 1Q Median 3Q Max
-26.3007 -5.9010 0.2208 5.8581 20.0910
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 86.94105 3.93131 22.11 <2e-16 ***
edat 0.75584 0.06857 11.02 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.241 on 98 degrees of freedom
Multiple R-squared: 0.5536, Adjusted R-squared: 0.549
F-statistic: 121.5 on 1 and 98 DF, p-value: < 2.2e-16Intercept ( _0) : pressió estimada si edat fos 0 (no sempre té sentit real)
Pendent (Estimate, _1) = : cada any d’edat incrementa la pressió mitjana en ~0.7 mmHg
Valor p: comprova si l’associació és significativa
R2 : proporció de variabilitat explicada pel model
En aquest cas concret (Intercept) = 86.94 és la pressió arterial estimada per una persona de 0 anys. No té sentit clínic, però és necessària per definir la recta.
Pendent (edat = 0.755)→ Per cada any d’edat, la pressió arterial augmenta en 0.71 mmHg, de mitjana. És el canvi mitjà en Y (pressió) per una unitat d’X (edat).
El valor p < 0.001 ens diu que l’associació entre edat i pressió és estadísticament significativa.
Per tal que els resultats del model (p-valor, intervals de confiança, etc.) siguin vàlids, cal que es compleixin els següents supòsits:
Linealitat: La relació entre XXX i YYY ha de ser lineal.
Normalitat dels residus: Els residus han de seguir una distribució normal. Ho podem veure amb un QQ-plot o un test de Shapiro.
Homocedasticitat (variància constant dels residus): Els residus han de tenir una variància constant al llarg dels valors de XXX.
Absència de punts influents: Cal detectar observacions que poden tenir una gran influència en el model. Comprova la distància de Cook:
No totes les observacions tenen el mateix impacte en l’ajust del model. Aquests punts s’anomenen valors influents i els detectem amb una mesura anomenada distància de Cook.
📌 Una observació amb Cook > 1 pot tenir una influència considerable. Però també val la pena revisar punts amb Cook més baixos però que destaquen visualment.
#comprobación grafica de linealidad
ggplot(dades, aes(x = edat, y = pressio)) +
geom_point(color = "gray40") +
geom_smooth(method = "lm", se = TRUE, color = "blue") +
labs(title = "Pressió arterial en funció de l'edat", x = "Edat (anys)", y = "Pressió sistòlica (mmHg)") `geom_smooth()` using formula = 'y ~ x'
residus <- residuals(model)
#Comprobacion de normalidad de residuus
## QQ-plot
qqnorm(residus); qqline(residus, col = "red") 
## Test de Shapiro-Wilk
shapiro.test(residus)
Shapiro-Wilk normality test
data: residus
W = 0.98779, p-value = 0.4932#Comprobacion de homocedasticidad de residus
plot(model$fitted.values, model$residuals,
xlab = "Valors ajustats", ylab = "Residus",
main = "Comprovació d'homocedasticitat")
abline(h = 0, col = "red")
#Comprobación de punts influents
plot(model, which = 4) # Cook's distance diagram
Amb mostres grans, el test de Shapiro pot detectar petites desviacions sense rellevància pràctica. El QQ-plot és sovint més útil visualment.
R genera quatre gràfics (per defecte) amb la funció plot(model) en aquest ordre:
1- Residus vs valors ajustats: comprova linealitat i homocedasticitat. 2- QQ-plot dels residus – comprova normalitat. 3- Scale-Location plot – variància constant dels residus. 4- Residuals vs leverage (Cook’s distance) – detecta punts influents.
| Supòsit | Comprovació visual | Test o eina |
|---|---|---|
| Linealitat | Diagrama de dispersió + línia lm |
— |
| Independència | Disseny de l’estudi | — |
| Normalitat dels residus | QQ-plot | shapiro.test() |
| Homocedasticitat | Residus vs valors ajustats | — |
| Absència d’influència | Cook’s distance | influence.measures() |
Ara que hem comprovat que el model compleix els supòsits bàsics, analitzem com de bé s’ajusta als nostres dades.
Es la proporció de la variabilitat de Y explicada per X.
summary(model)$r.squared # R2[1] 0.5535698summary(model)$adj.r.squared # R2 ajustat[1] 0.5490144EN el nostre model d’exemple, el 55.3% de la variabilitat de la pressió arterial s’explica per l’edat. El R2 ajustat corregeix el coeficient pel número de predictors.
L’AIC (Akaike Information Criterion) i el BIC (Bayesian Information Criterion) són criteris per valorar models, penalitzant la complexitat.
AIC(model)[1] 732.4929BIC(model)[1] 740.3084Com més baix sigui l’AIC o BIC, millor és el model. Només tenen significat relatiu: s’utilitzen per comparar diversos models ajustats a les mateixes dades. BIC penalitza més la complexitat que AIC.
📌 En aquest capítol, com que només tenim un model, no els interpretarem en profunditat, però són molt útils quan fem regressió múltiple o selecció de variables.
Aquest gràfic permet veure com de bones són les prediccions del model. Si tot fos perfecte, els punts estarien exactament sobre la línia vermella.
ggplot(data = dades, aes(x = fitted(model), y = pressio)) +
geom_point(color = "gray40") +
geom_abline(intercept = 0, slope = 1, col = "red", linetype = "dashed") +
labs(title = "Valors observats vs valors ajustats",
x = "Prediccions del model", y = "Pressió observada")
Quan tenim més d’un predictor que pot influir sobre una variable resposta quantitativa, utilitzem la regressió lineal múltiple.
Aquest model permet estimar l’efecte de cada variable predictora sobre la resposta, ajustant per les altres. És molt útil quan hi ha possibles confusors o quan volem fer una predicció més precisa.
Objectius del model múltiple
\[ Y = \beta\_0 + \beta\_1X_1 + \beta\_2X_2 + \dots + \beta\_pX_p + \varepsilon \] On:
📌 Aquest model és la base per a moltes tècniques avançades (regressió logística, models lineals generalitzats, etc.). El proper pas serà ajustar-lo amb lm() i interpretar-lo amb detall.
Suposem que volem predir la pressió arterial sistòlica d’un pacient utilitzant l’edat, el pes i el consum de tabac (fuma).
A diferència del model simple, ara estimarem l’efecte de cada variable sobre la pressió arterial, controlant per les altres.
# generem les dades de prova
set.seed(123)
n <- 100
edat <- sample(30:80, n, replace = TRUE)
pes <- rnorm(n, mean = 70, sd = 12)
fuma <- sample(c(0,1), n, replace = TRUE)
pressio <- 90 + 0.4 * edat + 0.3 * pes + 5 * fuma + rnorm(n, 0, 10)
dades2 <- data.frame(pressio, edat, pes, fuma = factor(fuma, labels = c("No", "Sí")))Ara que tenim més d’un predictor (edat, pes i si fuma), ajustem el model amb la funció lm() utilitzant dades2:
model_mult <- lm(pressio ~ edat + pes + fuma, data = dades2)
summary(model_mult)
Call:
lm(formula = pressio ~ edat + pes + fuma, data = dades2)
Residuals:
Min 1Q Median 3Q Max
-33.278 -7.769 1.992 7.352 28.922
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.33731 8.11307 11.874 < 2e-16 ***
edat 0.40038 0.08162 4.905 3.8e-06 ***
pes 0.21278 0.09973 2.134 0.0354 *
fumaSí 3.56173 2.22931 1.598 0.1134
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.95 on 96 degrees of freedom
Multiple R-squared: 0.2623, Adjusted R-squared: 0.2392
F-statistic: 11.38 on 3 and 96 DF, p-value: 1.89e-06Suposem que obtenim un resum com aquest:
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 91.215 5.124 17.8 <2e-16 ***
edat 0.37 0.07 5.3 <0.001 ***
pes 0.31 0.09 3.4 0.0011 **
fumaSí 5.02 2.15 2.3 0.023 *Intercept: pressió estimada d’una persona que no fuma, amb pes = 0 i edat = 0. No té interpretació clínica directa, però és necessària per definir el model.
edat: cada any d’edat incrementa la pressió en 0.37 mmHg, ajustant per pes i tabac.
pes: cada kg afegeix 0.31 mmHg, ajustant per edat i tabac.
fumaSí: les persones que fumen tenen, de mitjana, 5 mmHg més de pressió que les que no fumen, ajustant per edat i pes.
📌 Aquesta interpretació és condicional: cada efecte és estimat ajustant per la resta de variables.
summary(model_mult)$r.squared[1] 0.2622832summary(model_mult)$adj.r.squared[1] 0.2392295ggplot(data = dades2, aes(x = fitted(model_mult), y = pressio)) +
geom_point(color = "gray40") +
geom_abline(intercept = 0, slope = 1, col = "red", linetype = "dashed") +
labs(title = "Valors observats vs ajustats (model múltiple)",
x = "Prediccions del model", y = "Pressió observada")
📌 Si els punts estan prop de la línia diagonal, el model ajusta bé.
Si vols comparar aquest model amb el de regressió simple (només edat):
AIC(model)[1] 732.4929AIC(model_mult)[1] 768.3199Un AIC més baix indica un model millor. També pots comparar ( R^2 ) o fer una anàlisi jeràrquica amb anova(model, model_mult).
Aquest model ens permet fer interpretacions més realistes i ajustar per factors confusors.
El proper pas serà validar-lo i explorar interaccions si cal.
Per tal que els resultats del model siguin vàlids, cal que es compleixin els supòsits clàssics de la regressió lineal. Són els mateixos que en regressió simple, però ara n’afegim un de nou: la col·linealitat.
| Supòsit | Què implica | Com es comprova |
|---|---|---|
| Linealitat | Relació lineal entre cada predictor i la resposta | Gràfics, residus vs ajustats |
| Normalitat dels residus | Els errors segueixen una distribució normal | QQ-plot, shapiro.test() |
| Homocedasticitat | Variància constant dels residus | Residus vs ajustats |
| Absència de punts influents | Cap observació domina l’ajust del model | Cook’s distance |
| Col·linealitat | Els predictors no estan excessivament correlacionats entre ells | car::vif() |
par(mfrow = c(2,2))
plot(model_mult)
Aquests quatre gràfics et permeten detectar:
shapiro.test(residuals(model_mult))
Shapiro-Wilk normality test
data: residuals(model_mult)
W = 0.98306, p-value = 0.2284Col·linealitat vol dir que dues o més variables predictors estan molt correlacionades entre elles. Això pot fer que:
📌 No afecta la capacitat predictiva del model, però sí la interpretabilitat.
Per mesurar-la utilitzem el VIF (Variance Inflation Factor):
library(car)
vif(model_mult) edat pes fuma
1.009844 1.006812 1.003182 plot(model_mult, which = 4) # Cook's distance
📌 Punts amb Cook’s distance elevat poden tenir un efecte desproporcionat en el model. No cal eliminar-los directament, però cal revisar-los.
Aquestes comprovacions et permeten confirmar que el model és fiable abans d’interpretar-lo o utilitzar-lo per predicció.
El proper pas pot ser afegir interaccions o explorar models alternatius.
Fins ara hem suposat que cada variable afecta la resposta de forma independent. Però sovint, l’efecte d’un predictor depèn del valor d’un altre: això s’anomena interacció.
L’efecte del tabac sobre la pressió arterial és igual en totes les edats?
Potser fumar augmenta més la pressió en persones grans que en joves. Aquesta idea la representem amb una interacció entre edat i fuma.
model_inter <- lm(pressio ~ edat * fuma + pes, data = dades2)
summary(model_inter)
Call:
lm(formula = pressio ~ edat * fuma + pes, data = dades2)
Residuals:
Min 1Q Median 3Q Max
-33.421 -7.915 1.872 7.136 28.310
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 97.42428 9.16380 10.631 < 2e-16 ***
edat 0.38434 0.10268 3.743 0.000311 ***
fumaSí 1.07496 9.83126 0.109 0.913162
pes 0.20988 0.10083 2.081 0.040078 *
edat:fumaSí 0.04436 0.17075 0.260 0.795597
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11 on 95 degrees of freedom
Multiple R-squared: 0.2628, Adjusted R-squared: 0.2318
F-statistic: 8.467 on 4 and 95 DF, p-value: 6.917e-06La notació edat * fuma inclou:
edatfumaedat:fumaSuposem que veiem això:
Coefficients:
(Intercept) = 90.1
edat = 0.35
fumaSí = 4.0
edat:fumaSí = 0.15Això vol dir que:
ggplot(dades2, aes(x = edat, y = pressio, color = fuma)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE) +
labs(title = "Interacció entre edat i consum de tabac",
x = "Edat", y = "Pressió arterial")
Aquest gràfic mostra dues línies:
📌 Si les línies tenen pendents diferents, hi ha interacció.
Les interaccions ens ajuden a detectar efectes modulats entre variables. Però també poden complicar la interpretació: cal revisar si són plausibles i útils abans d’incloure-les sempre.
Quan construïm un model de regressió múltiple, sovint tenim moltes variables candidates. Però incloure-les totes pot fer el model:
Per això, busquem un model parsimoniós: amb poques variables, però útil i estable.
Hi ha diversos enfocaments. Els més habituals són:
step()Apliquem la selecció pas a pas basada en AIC:
model_tot <- lm(pressio ~ edat + pes + fuma + edat:fuma, data = dades2)
model_base <- lm(pressio ~ 1, data = dades2)
model_step <- step(model_base,
scope = list(lower = model_base, upper = model_tot),
direction = "both")Start: AIC=506.95
pressio ~ 1
Df Sum of Sq RSS AIC
+ edat 1 3233.3 12361 485.72
+ pes 1 787.2 14808 503.77
+ fuma 1 432.6 15162 506.14
<none> 15595 506.95
Step: AIC=485.72
pressio ~ edat
Df Sum of Sq RSS AIC
+ pes 1 551.0 11810 483.16
+ fuma 1 311.4 12050 485.17
<none> 12361 485.72
- edat 1 3233.3 15595 506.95
Step: AIC=483.16
pressio ~ edat + pes
Df Sum of Sq RSS AIC
+ fuma 1 305.90 11504 482.53
<none> 11810 483.16
- pes 1 551.04 12361 485.72
- edat 1 2997.12 14808 503.77
Step: AIC=482.53
pressio ~ edat + pes + fuma
Df Sum of Sq RSS AIC
<none> 11504 482.53
- fuma 1 305.90 11810 483.16
+ edat:fuma 1 8.17 11496 484.46
- pes 1 545.57 12050 485.17
- edat 1 2883.39 14388 502.90lower: el model mínim (sense predictors)upper: el model màxim (totes les variables)direction = "both": fa stepwise (pot afegir i treure)📌 Pots canviar a "forward" o "backward" si vols restringir el mètode.
Un cop seleccionat el model final, torna a validar els supòsits i interpreta els coeficients amb cura.
Un cop hem ajustat i seleccionat un model, cal validar-lo per assegurar-nos que:
Es divideixen les dades en dues parts:
set.seed(123)
n <- nrow(dades2)
idx <- sample(1:n, size = 0.7 * n)
training <- dades2[idx, ]
test <- dades2[-idx, ]
model_train <- lm(pressio ~ edat + pes + fuma, data = training)
pred_test <- predict(model_train, newdata = test)
# Error de predicció (RMSE)
rmse <- sqrt(mean((test$pressio - pred_test)^2))
rmse[1] 9.425799📌 RMSE (Root Mean Squared Error): quantifica com s’equivoca el model. Com més petit, millor.
Fa diverses divisions de training/test i calcula l’error mitjà.
library(boot)
model_cv <- glm(pressio ~ edat + pes + fuma, data = dades2)
cv.err <- cv.glm(dades2, model_cv, K = 10)
cv.err$delta[1] 123.0042 122.5808delta[1]: error mitjà d’entrenamentdelta[2]: error estimat d’una predicció nova (més fiable)📌 La cross-validation és més estable i precisa que una sola partició.
Alternativa per validar la variabilitat dels coeficients:
library(boot)
boot_fn <- function(data, i) {
coef(lm(pressio ~ edat + pes + fuma, data = data[i, ]))
}
resultats_boot <- boot(data = dades2, statistic = boot_fn, R = 1000)
boot.ci(resultats_boot, type = "bca", index = 2) # interval bootstrap per edatBOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates
CALL :
boot.ci(boot.out = resultats_boot, type = "bca", index = 2)
Intervals :
Level BCa
95% ( 0.2416, 0.5710 )
Calculations and Intervals on Original Scale📌 El bootstrap ajuda a veure com canvien els coeficients si repetim el mostreig. Ens dona intervals de confiança més robustos.
| Mètode | Objectiu principal |
|---|---|
| Training/Test | Predicció a noves dades (ràpid i visual) |
| Cross-validation | Estimació precisa d’error de predicció |
| Bootstrap | Fiabilitat dels coeficients |
Tots tres mètodes són útils i complementaris.
Un bon model hauria de preveure bé i tenir coeficients estables.
Un cop hem ajustat, interpretat, validat i revisat el model, podem fer una síntesi estructurada dels resultats.
Estimar la pressió arterial sistòlica a partir de tres variables predictives:
El model final inclou:
summary(model_mult)$coefficients Estimate Std. Error t value Pr(>|t|)
(Intercept) 96.3373129 8.11306741 11.874339 1.515568e-20
edat 0.4003836 0.08162492 4.905164 3.804656e-06
pes 0.2127800 0.09972502 2.133667 3.541881e-02
fumaSí 3.5617252 2.22930631 1.597683 1.134000e-01Totes les variables són estadísticament significatives i clínicament rellevants.
summary(model_mult)$r.squared[1] 0.2622832summary(model_mult)$adj.r.squared[1] 0.2392295El model explica aproximadament el 55% de la variabilitat de la pressió arterial.
Aquest valor és acceptable en dades reals, però indica que hi ha altres factors implicats.
rmse[1] 9.425799cv.err$delta[2][1] 122.5808L’error de predicció és moderadament baix, el que indica que el model generalitza raonablement bé.
L’anàlisi de la variància (ANOVA) és una tècnica estadística per comparar la mitjana d’una variable quantitativa entre dos o més grups.
Tot i que sovint es presenta com una tècnica separada, ANOVA és en realitat un cas particular de regressió lineal, on el predictor és categòric.
Analitzarem si hi havia diferències d’edat mitjana entre passatgers de Primera, Segona i Tercera classe.
# Seleccionem variables d'interès i netegem NA
library(readxl)
library(dplyr)
library(ggplot2)
dades_titanic <- titanic %>%
select(Survived, Pclass, Sex, Age) %>%
filter(!is.na(Age))
dades_titanic$Pclass <- factor(dades_titanic$Pclass,
levels = c(1,2,3),
labels = c("Primera", "Segona", "Tercera"))
### Visualització prèvia
ggplot(dades_titanic, aes(x = Pclass, y = Age)) +
geom_boxplot(fill = "lightblue") +
labs(title = "Edat segons classe del passatger",
x = "Classe", y = "Edat")
Aquest gràfic suggereix que hi ha diferències clares en l’edat entre les classes. Comprovarem si les petites diferències visualitzades, són estadísticament significatives.
model_anova <- lm(Age ~ Pclass, data = dades_titanic)
anova(model_anova)Analysis of Variance Table
Response: Age
Df Sum Sq Mean Sq F value Pr(>F)
Pclass 2 20930 10464.8 57.444 < 2.2e-16 ***
Residuals 711 129527 182.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1La sortida ens mostra el contrast F que indica si hi ha diferències significatives entre les mitjanes dels grups. Igual que en regressió lineal, cal que es compleixin aquests supòsits:
| Supòsit | Comprovació | Codi |
|---|---|---|
| Linealitat | Ja garantida per variable categòrica | — |
| Normalitat dels residus | QQ-plot i test de Shapiro | shapiro.test() |
| Homocedasticitat | Variància similar entre grups | bartlett.test() o leveneTest() |
# Normalitat dels residus
shapiro.test(residuals(model_anova))
Shapiro-Wilk normality test
data: residuals(model_anova)
W = 0.99146, p-value = 0.0003856# Homocedasticitat
bartlett.test(Age ~ Pclass, data = dades_titanic)
Bartlett test of homogeneity of variances
data: Age by Pclass
Bartlett's K-squared = 7.818, df = 2, p-value = 0.02006##
plot(model_anova)
📌 No obstant al QQ-plot, s’observa una bona normalitat essent els extrems i allunyats de la línia teòrica.
En el nostre exemple, el test de Shapiro-Wilk per comprovar la normalitat dels residus ha donat un resultat significatiu, però el QQ-plot mostra una forma força acceptable.
Això passa sovint quan tenim moltes observacions: el test és molt sensible i pot detectar desviacions petites que no són rellevants.
“Si el QQ-plot es veu raonablement bé i el model té molts casos (n > 50), pots continuar amb el model però indicant que hi ha una possible desviació lleu.”
NO obstant si sospites o clarament veus que el model ANOVA no compleix els supòsits i per tant:
Aquest test no assumeix normalitat ni homocedasticitat. Comprova si almenys un grup té una mediana diferent.
kruskal.test(Age ~ Pclass, data = dades_titanic)
Kruskal-Wallis rank sum test
data: Age by Pclass
Kruskal-Wallis chi-squared = 95.995, df = 2, p-value < 2.2e-16En aquest cas el test no paramètric també surt molt significatiu i per tant podem asegurar que hi han diferències en edat per classe de pasatge.
summary(model_anova)
Call:
lm(formula = Age ~ Pclass, data = dades_titanic)
Residuals:
Min 1Q Median 3Q Max
-37.313 -7.878 -0.878 7.859 48.859
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 38.2334 0.9897 38.633 < 2e-16 ***
PclassSegona -8.3558 1.4257 -5.861 7.04e-09 ***
PclassTercera -13.0928 1.2217 -10.717 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.5 on 711 degrees of freedom
Multiple R-squared: 0.1391, Adjusted R-squared: 0.1367
F-statistic: 57.44 on 2 and 711 DF, p-value: < 2.2e-16Els coeficients indiquen les diferències d’edat mitjana respecte al grup de referència (Primera classe):
Si el contrast F és significatiu, podem fer comparacions entre grups amb correcció per múltiples comparacions:
TukeyHSD(aov(model_anova)) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = model_anova)
$Pclass
diff lwr upr p adj
Segona-Primera -8.355811 -11.704133 -5.007489 0.0000000
Tercera-Primera -13.092821 -15.962198 -10.223445 0.0000000
Tercera-Segona -4.737010 -7.676279 -1.797742 0.0004884Això ens indica quins parells de grups difereixen significativament.
Per fer comparacions múltiples no paramètriques, es fa servir el test de Dunn. FSA::dunnTest()
library(FSA)Registered S3 methods overwritten by 'FSA':
method from
confint.boot car
hist.boot car ## FSA v0.10.0. See citation('FSA') if used in publication.
## Run fishR() for related website and fishR('IFAR') for related book.
S'està adjuntant el paquet: 'FSA'L'objecte següent està emmascarat per 'package:car':
bootCasedunnTest(Age ~ Pclass, data = dades_titanic, method = "holm")Dunn (1964) Kruskal-Wallis multiple comparison p-values adjusted with the Holm method. Comparison Z P.unadj P.adj
1 Primera - Segona 4.798551 1.598180e-06 3.196359e-06
2 Primera - Tercera 9.762001 1.638887e-22 4.916662e-22
3 Segona - Tercera 4.063515 4.833916e-05 4.833916e-05Quan la variable dependent és binària (sí/no, viu/mor, positiu/negatiu), la regressió lineal no és adequada, ja que pot donar valors fora de l’interval [0,1].
La regressió logística ens permet modelitzar la probabilitat d’un esdeveniment, ajustant-la a partir de variables explicatives quantitatives o categòriques.
| Model | Variable resposta | Resultat estimat | Exemple |
|---|---|---|---|
| Regressió lineal | Quantitativa | Valor mitjà | Pressió arterial, colesterol |
| Regressió logística | Binària (0/1) | Probabilitat (entre 0 i 1) | Supervivència, diagnòstic, resposta a tractament |
En regressió logística, el model ajusta el logit de la probabilitat:
[ (p) = () = _0 + _1 X ]
Això vol dir que la relació entre X i la probabilitat no és lineal, sinó en forma de S (sigmoide).
A partir d’aquí, ajustarem un model amb glm() per predir la probabilitat de supervivència (Survived) al Titanic a partir de variables com sexe, edat i classe.
Utilitzarem la funció glm() amb la família binomial per modelitzar la probabilitat de supervivència (Survived) segons:
titanic <- read.csv("titanic.csv")
library(dplyr)
dades_titanic <- titanic %>%
select(Survived, Pclass, Sex, Age) %>%
filter(!is.na(Age)) %>%
mutate(
Survived = factor(Survived, levels = c(0,1), labels = c("No", "Sí")),
Pclass = factor(Pclass, levels = c(1,2,3), labels = c("Primera", "Segona", "Tercera")),
Sex = factor(Sex)
)model_log <- glm(Survived ~ Sex + Pclass + Age,
data = dades_titanic, family = binomial)
summary(model_log)
Call:
glm(formula = Survived ~ Sex + Pclass + Age, family = binomial,
data = dades_titanic)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.777013 0.401123 9.416 < 2e-16 ***
Sexmale -2.522781 0.207391 -12.164 < 2e-16 ***
PclassSegona -1.309799 0.278066 -4.710 2.47e-06 ***
PclassTercera -2.580625 0.281442 -9.169 < 2e-16 ***
Age -0.036985 0.007656 -4.831 1.36e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 964.52 on 713 degrees of freedom
Residual deviance: 647.28 on 709 degrees of freedom
AIC: 657.28
Number of Fisher Scoring iterations: 5Els coeficients del model són en escala log-odds. Per interpretar-los com a odds ratio, els transformem amb exp().
En regressió logística, el model s’ajusta sobre l’escala del logaritme dels odds:
[ (p) = () = _0 + _1X_1 + + _kX_k ]
Això vol dir que els coeficients del model (summary(model_log)) representen:
el canvi en el logaritme dels odds (log-odds) per cada unitat de canvi en el predictor, ajustant per la resta.
Però com que els log-odds no són fàcils d’interpretar, el que fem habitualment és exponenciar els coeficients:
exp(coef(model_log)) (Intercept) Sexmale PclassSegona PclassTercera Age
43.68534331 0.08023617 0.26987422 0.07572664 0.96369033 Això ens dona les odds ratio (OR), que sí que són interpretables de forma més directa:
Si
exp(β₁) = 2, això vol dir que els odds de sobreviure es dupliquen per cada unitat de canvi en la variable ( X_1 ), mantenint la resta constant.
Per a variables categòriques, la interpretació és per comparació amb el grup de referència.
exp(confint(model_log)) # Intervals de confiança de les odds ratioS'està esperant que es faci l'anàlisi de rendiment... 2.5 % 97.5 %
(Intercept) 20.37890724 98.3863313
Sexmale 0.05293643 0.1194848
PclassSegona 0.15515074 0.4621731
PclassTercera 0.04299250 0.1297997
Age 0.94905346 0.9780124Aquest enfocament ens permet fer una interpretació clínica clara dels efectes de cada variable.
exp(coef(model_log)) # Odds Ratio (Intercept) Sexmale PclassSegona PclassTercera Age
43.68534331 0.08023617 0.26987422 0.07572664 0.96369033 exp(confint(model_log)) # Interval de confiança dels ORS'està esperant que es faci l'anàlisi de rendiment... 2.5 % 97.5 %
(Intercept) 20.37890724 98.3863313
Sexmale 0.05293643 0.1194848
PclassSegona 0.15515074 0.4621731
PclassTercera 0.04299250 0.1297997
Age 0.94905346 0.9780124L’odds ratio per
Sexe = donaés 4.5
→ Això vol dir que les dones tenen 4.5 vegades més probabilitat d’haver sobreviscut que els homes, ajustant per edat i classe.
Tot i que no cal comprovar normalitat ni homocedasticitat com en la regressió lineal, la regressió logística té altres supòsits importants:
| Supòsit | Explicació breu | Com es comprova |
|---|---|---|
| 1. Variable dependent binària | La resposta ha de ser 0/1 (o factor amb dos nivells) | Ja ho hem fet amb Survived |
| 2. Observacions independents | No hi ha dades repetides ni agrupades | Revisió del disseny |
| 3. Linealitat del logit | Relació lineal entre els predictors quantitatius i el logit(p) | Test de Box-Tidwell o visualment |
| 4. Absència de col·linealitat | Els predictors no han d’estar fortament correlacionats | car::vif() |
| 5. Mostra prou gran | Per estabilitat i validesa dels coeficients | n ≥ 10 esdeveniments per predictor |
Només afecta als predictors quantitatius (com Age en el nostre cas).
Una manera senzilla d’inspeccionar-ho és representar:
dades_titanic$prob_pred <- predict(model_log, type = "response")
dades_titanic <- dades_titanic %>%
mutate(logit = log(prob_pred / (1 - prob_pred)))
ggplot(dades_titanic, aes(x = Age, y = logit)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "loess", se = FALSE, col = "red") +
labs(title = "Linealitat del logit respecte a l'edat")`geom_smooth()` using formula = 'y ~ x'
📌 Si la línia vermella fa una corba, potser caldria transformar Age o categoritzar-la.
{r} library(car) vif(model_log)}`
Aquests supòsits no són difícils de revisar, però són essencials per garantir la validesa del model.
Per validar un model de regressió logística, ens fixem sobretot en la seva capacitat de discriminació, és a dir:
📌 Com bé és capaç el model de distingir entre els que sobreviuen i els que no.
Un cop ajustat el model, podem obtenir les probabilitats predites de supervivència:
dades_titanic$prob_pred <- predict(model_log, type = "response")
head(dades_titanic$prob_pred)[1] 0.10526285 0.91463372 0.55842450 0.92290788 0.06780678 0.32235446Aquest valor és la probabilitat estimada de Survived = Sí per a cada individu.
La corba ROC representa:
Com més s’apropa la corba a la cantonada superior esquerra, millor.
library(pROC)
roc_obj <- roc(dades_titanic$Survived, dades_titanic$prob_pred)
plot(roc_obj, col = "blue", main = "Corba ROC - Model logístic")
auc(roc_obj)Area under the curve: 0.8523L’AUC (Area Under the Curve) ens dona una única xifra per quantificar la discriminació:
| AUC | Interpretació |
|---|---|
| 0.5 | Aleatori (el model no discrimina) |
| 0.6–0.7 | Discriminació feble |
| 0.7–0.8 | Bona discriminació |
| 0.8–0.9 | Molt bona discriminació |
| > 0.9 | Excel·lent (pot ser sobreajustat!) |
Exemple: un AUC de 0.83 vol dir que el model té un 83% de probabilitat d’ordenar correctament un cas positiu per damunt d’un negatiu.
També es pot:
Amb aquesta validació, podem afirmar si el model és útil per predir esdeveniments i prendre decisions.
Al llarg d’aquest capítol hem après que un model estadístic és una eina per descriure, predir o explicar una variable resposta a partir d’unes variables explicatives.
Hem vist que molts dels models més habituals formen part d’una mateixa família: el model lineal.
Depèn de:
| Model | Variable dependent | Predictors | Exemple |
|---|---|---|---|
| Regressió lineal simple | Quantitativa | 1 variable quantitativa | Pressió ~ edat |
| Regressió lineal múltiple | Quantitativa | 2+ variables (qualsevol tipus) | Pressió ~ edat + pes + sexe |
| ANOVA d’un factor | Quantitativa | 1 variable categòrica | Edat ~ classe |
| Regressió logística | Binària (Sí/No) | 1 o més predictors | Supervivència ~ sexe + classe + edat |
Si vols extreure resultats de manera ordenada i convertir-los en taules fàcilment exportables, pots utilitzar la funció broom::tidy(model). És especialment útil per informes o quan treballes amb molts models.
Tanmateix, per comprendre bé el funcionament d’un model, és millor fer primer la interpretació pas a pas amb summary(), coef(), exp() i gràfics.