Curs d’Estadística amb R per a Metges

Capítol 4: Models estadístics en biomedicina

Autor/a

P. Loma-Osorio

Publicat

15 Octubre, 2026

1 Capítol 4: Models estadístics en biomedicina

Codi 
library(tidyverse)
library(car)
library(boot)
library(pROC)
library(gtsummary)

1.1 Què és un model?

Un model estadístic és una representació matemàtica simplificada de la realitat que descriu com es relaciona una variable clínica de resposta (variable dependent, \(Y\)) amb un o diversos factors d’estudi o confusors (variables independents o predictors, \(X\)):

\[Y = f(X) + \varepsilon\]

On: * \(Y\): variable dependent o resposta clínica (ex: pressió arterial, mortalitat). * \(X\): variable o variables independents (predictors o factors de risc). * \(\varepsilon\): error aleatori o soroll biològic no explicat.

En l’àmbit biomèdic utilitzem els models amb tres finalitats essencials: 1. Explicar: Quantificar l’associació directa entre factors i desenllaços clínics. 2. Controlar confusió: Aïllar l’efecte d’una variable ajustant per la resta. 3. Predir: Estimar valors o probabilitats futures per ajudar en la presa de decisions.

1.1.1 La mitjana com a model: introducció intuïtiva a la regressió

Imaginem que volem estimar l’alçada d’un pacient que entra per la porta de consultes, però no disposem de cap informació clínica addicional.

En absència de qualsevol dada, el millor model possible és la mitjana mostral (\(\bar{Y}\)):

\[Y_i = \bar{Y} + \varepsilon_i\]

Aquest model elemental diu: “Com que no sé res més d’aquest subjecte, estimo que mesura exactament la mitjana del grup”.

Codi 
set.seed(42)
alçades <- rnorm(100, mean = 170, sd = 10)
mitjana <- mean(alçades)

plot(alçades, rep(1, 100), pch = 16, col = "gray50", yaxt = "n",
     xlab = "Alçada (cm)", ylab = "", 
     main = "La mitjana com a model (estimació basal constant)")
abline(v = mitjana, col = "firebrick", lwd = 2)
text(mitjana + 1, 1.2, paste("Mitjana =", round(mitjana, 1), "cm"), 
     col = "firebrick", pos = 4)
Figura 1

1.1.2 Els residus com a error del model

Com és evident, pràcticament cap pacient mesurarà exactament la mitjana. La diferència individual entre el valor observat real de cada pacient (\(Y_i\)) i l’estimació del model (\(\hat{Y}_i = \bar{Y}\)) s’anomena residu o error:

\[\varepsilon_i = Y_i - \hat{Y}_i\]

El conjunt dels residus ens indica com d’allunyat està el món real del nostre model:

Codi 
residus <- alçades - mitjana

hist(residus, col = "#BBDEFB", border = "white",
     main = "Distribució dels residus respecte a la mitjana",
     xlab = "Alçada real - Mitjana (cm)", 
     breaks = 10)
abline(v = 0, col = "firebrick", lwd = 2, lty = 2)
Figura 2

La desviació estàndard (SD) és una mesura de la magnitud mitjana d’aquests errors. Com més gran sigui la dispersió dels residus, menor és la precisió del nostre model.

1.1.3 Regressió lineal com a pas següent

Ara imaginem que sabem una dada addicional de cada pacient: la seva edat. Ja no ens cal assumir una mitjana plana horitzontal constant, sinó que podem ajustar una recta amb pendent que s’adapti als canvis d’edat:

Codi 
set.seed(123)
edat <- sample(20:80, 100, replace = TRUE)
alçada <- 150 + 0.4 * edat + rnorm(100, mean = 0, sd = 5)
dades <- data.frame(edat, alçada)

model_intro <- lm(alçada ~ edat, data = dades)

plot(dades$edat, dades$alçada, pch = 16, col = "gray40",
     xlab = "Edat (anys)", ylab = "Alçada (cm)",
     main = "Model de regressió lineal: adaptació segons predictor")
abline(h = mean(dades$alçada), col = "gray70", lty = 3, lwd = 1.5)
abline(model_intro, col = "#1976D2", lwd = 2.5)
legend("topleft", legend = c("Model mitjana", "Model lineal"),
       col = c("gray70", "#1976D2"), lty = c(3, 1), lwd = c(1.5, 2.5), bty = "n")
Figura 3

En analitzar la variabilitat tenim tres components: 1. Variabilitat total: Quant varien les respostes observades al voltant de la mitjana. 2. Variabilitat explicada: La fracció que aconsegueix capturar la recta de regressió a través del predictor. 3. Variabilitat no explicada: Els residus individuals que encara queden al voltant de la línia.

1.2 Regressió lineal simple

En la regressió lineal simple modelitzem la relació entre dues variables quantitatives:

\[Y = \beta_0 + \beta_1 X + \varepsilon\]

On: * \(Y\): variable dependent (resposta contínua). * \(X\): variable independent (predictora quantitativa). * \(\beta_0\): interceptació (valor esperat de \(Y\) quan \(X = 0\)). * \(\beta_1\): pendent o slope (canvi esperat en \(Y\) per cada unitat d’increment en \(X\)). * \(\varepsilon\): error aleatori o residu.

1.2.1 Ajust del model

Analitzem la relació entre l’edat i la pressió arterial sistòlica (TAS) en una cohort simulada de 100 pacients:

Codi 
set.seed(123)
n <- 100
edat <- sample(30:80, n, replace = TRUE)
pressio <- 90 + 0.7 * edat + rnorm(n, mean = 0, sd = 10)
dades <- data.frame(edat, pressio)

# Ajust del model de regressió lineal simple
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-16

1.2.2 Interpretació dels coeficients

  • Intercept (\(\beta_0 = 86.94\) mmHg): Pressió arterial estimada per a una edat teòrica de 0 anys. No té significat biològic aïllat, però fixa l’origen matemàtic de la recta.
  • Pendent (\(\beta_1 = 0.756\) mmHg/any): Per cada any adicional d’edat del subjecte, la pressió arterial sistòlica mitjana s’incrementa en 0.76 mmHg.
  • Significació estadística: El contrast d’hipòtesi per al pendent (\(t = 11.02, p < 0.001\)) confirma que l’edat s’associa de forma estadísticament significativa a la pressió sistòlica.

1.2.3 Supòsits del model de regressió lineal simple

Perquè les conclusions del model siguin vàlides cal avaluar cinc supòsits essencials:

Supòsit Què avalua Com es comprova
Linealitat Relació lineal directa entre \(X\) i \(Y\) Diagrama de dispersió i residus vs. ajustats
Independència Observacions no aparellades ni agrupades Disseny de l’estudi
Normalitat dels residus Errors distribuïts segons una campana de Gauss QQ-plot i shapiro.test()
Homocedasticitat Variància constant dels residus al llarg dels valors Residus vs. ajustats i Scale-Location
Absència de punts influents Cap cas aïllat condiciona la recta Distància de Cook (plot(model, which = 4))

1.2.3.1 Comprovació gràfica de linealitat

Codi 
ggplot(dades, aes(x = edat, y = pressio)) + 
  geom_point(color = "gray40", alpha = 0.7) + 
  geom_smooth(method = "lm", se = TRUE, color = "#1976D2") +
  labs(title = "Pressió arterial en funció de l'edat", 
       x = "Edat (anys)", 
       y = "Pressió sistòlica (mmHg)") +
  theme_classic()
Figura 4

1.2.3.2 Detecció de punts influents (Distància de Cook)

Una observació amb una distància de Cook \(> 1\) pot exercir una influència desproporcionada sobre l’estimació dels paràmetres:

Codi 
plot(model, which = 4, col = "#D32F2F")
Figura 5

1.2.3.3 Els quatre gràfics diagnòstics estàndard d’R

Executant plot(model) obtenim la bateria diagnòstica completa:

Codi 
par(mfrow = c(2, 2))
plot(model)
par(mfrow = c(1, 1))
Figura 6
  1. Residuals vs Fitted: Avalua linealitat i homocedasticitat (la línia vermella ha de ser plana al voltant de 0).
  2. Normal Q-Q: Avalua normalitat dels residus (els punts s’han de situar sobre la diagonal).
  3. Scale-Location: Avalua homocedasticitat amb residus estandarditzats.
  4. Residuals vs Leverage: Detecta punts atípics i palanquejament (leverage).

1.2.4 Diagnòstic del model: Bondat d’ajust

1.2.4.1 Coeficient de determinació (\(R^2\))

Quantifica la fracció de variància de la variable dependent que queda explicada pel model:

Codi 
summary(model)$r.squared      # R2
[1] 0.5535698
Codi 
summary(model)$adj.r.squared  # R2 ajustat
[1] 0.5490144

En aquest model, el 55.4% de la variabilitat de la pressió arterial s’explica directament per l’edat del pacient.

1.2.4.2 Criteris d’informació: AIC i BIC

L’AIC (Akaike Information Criterion) i el BIC (Bayesian Information Criterion) penalitzen la complexitat per evitar el sobreajust. Valors inferiors indiquen millor equilibri entre parsimònia i precisió:

Codi 
AIC(model)
[1] 732.4929
Codi 
BIC(model)
[1] 740.3084

1.2.4.3 Gràfic de valors observats vs. ajustats

Permet confrontar la concordança directa entre les prediccions del model i les observacions reals:

Codi 
ggplot(data = dades, aes(x = fitted(model), y = pressio)) +
  geom_point(color = "gray40", alpha = 0.7) +
  geom_abline(intercept = 0, slope = 1, col = "red", linetype = "dashed", linewidth = 1) +
  labs(title = "Valors observats vs. ajustats",
       x = "Prediccions del model (mmHg)", 
       y = "Pressió observada (mmHg)") +
  theme_classic()
Figura 7

1.3 Regressió lineal múltiple

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 l’eina bàsica en biomedicina per avaluar associacions clíniques controlant per possibles factors de confusió o per construir sistemes de predicció.

Objectius del model múltiple: * Estimar l’efecte individual de cada predictor. * Controlar variables confusores (edat, comorbiditats, gravetat basal). * Fer prediccions ajustades a perfils concrets de pacients. * Comparar la contribució relativa de cada factor.

1.3.1 Estructura general del model

\[Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_p X_p + \varepsilon\]

On: * \(Y\): variable dependent o desenllaç continu. * \(X_1, X_2, \dots, X_p\): variables independents o predictors (continus o categòrics). * \(\beta_0\): interceptació (valor esperat de \(Y\) quan tots els predictors valen 0). * \(\beta_i\): efecte estimat del predictor \(X_i\) mantenint constants tots els altres predictors. * \(\varepsilon\): error aleatori o residu.

Suposem que volem predir la pressió arterial sistòlica d’un pacient utilitzant l’edat, el pes i el consum de tabac (fuma).

Codi 
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, mean = 0, sd = 10)

dades2 <- data.frame(
  pressio = pressio, 
  edat = edat, 
  pes = pes, 
  fuma = factor(fuma, levels = c(0, 1), labels = c("No", "Sí"))
)

1.3.2 Ajust del model de regressió múltiple

Ajustem el model amb la funció lm() incloent els tres predictors sumats:

Codi 
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-06

1.3.3 Interpretació dels coeficients

La interpretació en regressió múltiple és estrictament condicional: cada coeficient quantifica l’efecte d’aquella variable mantenint constants la resta de variables del model:

  • Intercept (\(\beta_0 = 96.34\) mmHg): Pressió arterial estimada per a una persona que no fuma, amb edat 0 i pes 0 kg. Fixa l’origen de coordenades del pla matemàtic.
  • edat (\(\beta_1 = 0.40\) mmHg/any; \(p < 0.001\)): Per cada any d’edat addicional, la pressió arterial s’incrementa de mitjana en 0.40 mmHg, després d’ajustar per pes i consum de tabac.
  • pes (\(\beta_2 = 0.21\) mmHg/kg; \(p = 0.035\)): Cada quilogram addicional de pes s’associa a un augment mitjà de 0.21 mmHg de pressió sistòlica, ajustant per edat i tabac.
  • fumaSí (\(\beta_3 = 3.56\) mmHg; \(p = 0.113\)): Els individus fumadors presenten, de mitjana, 3.56 mmHg més de pressió arterial sistòlica que els no fumadors després de controlar per edat i pes. En aquesta mostra de 100 subjectes, aquesta diferència no assoleix la significació estadística al llindar estàndard de 0.05.

1.3.4 Diagnòstics del model múltiple

1.3.4.1 Coeficient de determinació (\(R^2\))

Codi 
summary(model_mult)$r.squared      # R2 brut
[1] 0.2622832
Codi 
summary(model_mult)$adj.r.squared  # R2 ajustat
[1] 0.2392295
  • El model explica el 26.2% (\(R^2 = 0.262\)) de la variabilitat de la pressió arterial.
  • L’\(R^2\) ajustat (\(0.239\)) corregeix el valor segons el nombre de variables incloses, penalitzant l’addició de predictors que no aportin informació rellevant.

1.3.4.2 Valors observats vs. ajustats

Codi 
ggplot(data = dades2, aes(x = fitted(model_mult), y = pressio)) +
  geom_point(color = "gray40", alpha = 0.7) +
  geom_abline(intercept = 0, slope = 1, col = "red", linetype = "dashed", linewidth = 1) +
  labs(title = "Valors observats vs. ajustats (model múltiple)",
       x = "Prediccions del model (mmHg)", 
       y = "Pressió observada (mmHg)") +
  theme_classic()
Figura 8

1.3.4.3 Comparació amb el model simple

Podem contrastar si la inclusió de pes i tabac aporta una millora substancial comparant criteris d’informació:

Codi 
AIC(model)       # Model simple (només edat)
[1] 732.4929
Codi 
AIC(model_mult)  # Model múltiple (edat, pes, tabac)
[1] 768.3199

1.3.5 Validació dels supòsits en regressió múltiple

Supòsit Què implica Com es comprova
Linealitat Relació lineal entre cada predictor i la resposta contínua Gràfic de residus vs. ajustats
Normalitat dels residus Errors distribuïts de forma gaussiana QQ-plot i shapiro.test()
Homocedasticitat Variància constant dels residus Gràfic de residus vs. valors ajustats
Absència de punts influents Cap observació té un pes excessiu sobre els coeficients Distància de Cook (plot(model, which = 4))
Absència de col·linealitat Els predictors no estan fortament correlacionats entre ells Factor d’inflació de la variància (car::vif())

1.3.5.1 Gràfics de diagnòstic complets

Codi 
par(mfrow = c(2, 2))
plot(model_mult)
par(mfrow = c(1, 1))
Figura 9

1.3.5.2 Normalitat dels residus

Codi 
shapiro.test(residuals(model_mult))

    Shapiro-Wilk normality test

data:  residuals(model_mult)
W = 0.98306, p-value = 0.2284

1.3.5.3 Col·linealitat: explicació i detecció

La col·linealitat apareix quan dos o més predictors estan fortament correlacionats entre ells, inflant els errors estàndard. La detectem amb el VIF (Variance Inflation Factor):

Codi 
library(car)
vif(model_mult)
    edat      pes     fuma 
1.009844 1.006812 1.003182
  • VIF \(\approx 1\): Absència total de correlació entre predictors.
  • VIF \(> 5\): Col·linealitat moderada.
  • VIF \(> 10\): Col·linealitat severa que requereix intervenció.

1.3.5.4 Punts influents

Codi 
plot(model_mult, which = 4, col = "#D32F2F")
Figura 10

1.3.6 Interaccions entre predictors

En biomedicina és freqüent que l’efecte d’un factor depengui de la presència o nivell d’un altre (modificació de l’efecte o interacció).

La sintaxi edat * fuma avalua l’efecte principal de l’edat, l’efecte principal del tabac i el terme combinat d’interacció (edat:fuma):

Codi 
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-06

1.3.6.1 Visualització de la interacció

Codi 
ggplot(dades2, aes(x = edat, y = pressio, color = fuma)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE) +
  scale_color_manual(values = c("No" = "#1976D2", "Sí" = "#D32F2F")) +
  labs(title = "Interacció entre edat i consum de tabac sobre la pressió",
       x = "Edat (anys)", 
       y = "Pressió arterial sistòlica (mmHg)",
       color = "Fuma?") +
  theme_classic()
Figura 11

1.3.7 Selecció de variables en regressió lineal múltiple

L’objectiu és trobar un model parsimoniós: el menor nombre de predictors capaços de retenir la màxima informació clínica.

Codi 
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",
  trace = FALSE
)

summary(model_step)

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-06

Criteri clínic: Els algorismes automàtics de selecció mai han de substituir el raonament fisiopatològic o la literatura prèvia. Un factor amb valor pronòstic reconegut o un confusor conegut s’ha de mantenir al model encara que el seu valor p no assoleixi la significació estadística a la mostra.

1.3.8 Validació interna del model continu

1.3.8.1 1. Divisió Training / Test

Codi 
set.seed(123)
n_obs <- nrow(dades2)
idx <- sample(1:n_obs, size = 0.7 * n_obs)

training <- dades2[idx, ]
test <- dades2[-idx, ]

model_train <- lm(pressio ~ edat + pes + fuma, data = training)
pred_test <- predict(model_train, newdata = test)

# Càlcul de l'error quadràtic mitjà (RMSE)
rmse <- sqrt(mean((test$pressio - pred_test)^2))
rmse
[1] 9.425799

1.3.8.2 2. Validació creuada (10-Fold Cross-Validation)

Codi 
library(boot)

model_glm <- glm(pressio ~ edat + pes + fuma, data = dades2)
cv_err <- cv.glm(dades2, model_glm, K = 10)

# Error estimat de predicció generalitzable (MSE)
cv_err$delta[2]
[1] 122.5808

1.3.8.3 3. Bootstrap

Codi 
boot_fn <- function(data, indices) {
  coef(lm(pressio ~ edat + pes + fuma, data = data[indices, ]))
}

set.seed(123)
resultats_boot <- boot(data = dades2, statistic = boot_fn, R = 1000)

# Interval de confiança robust pel coeficient d'edat (índex 2)
boot.ci(resultats_boot, type = "bca", index = 2)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = resultats_boot, type = "bca", index = 2)

Intervals : 
Level       BCa          
95%   ( 0.2413,  0.5664 )  
Calculations and Intervals on Original Scale

1.3.9 Informe final del model de regressió múltiple

  • Objectiu: Identificar els determinants de la pressió arterial sistòlica i quantificar el seu pes independent.
  • Variables incloses: Edat (anys), pes corporal (kg) i hàbit tabàquic (Sí/No).
  • Bondat d’ajust: El model explica el 26.2% de la variabilitat de la pressió (\(R^2 = 0.262\), \(R^2 \text{ ajustat} = 0.239\)).
  • Validació interna: S’obté un error mitjà de predicció sobre la mostra de test de 9.43 mmHg (RMSE) i una estimació per validació creuada de 10 particions que confirma l’estabilitat de l’ajust.
  • Conclusions: L’edat i el pes s’associen de forma lineal, independent i estadísticament significativa amb la pressió arterial. Es verifiquen satisfactòriament els supòsits de linealitat, normalitat dels residus, homocedasticitat i absència de col·linealitat (tots els valors de VIF \(\approx 1\)).

1.4 ANOVA d’un factor com a model lineal

L’anàlisi de la variància (ANOVA) s’utilitza habitualment per comparar la mitjana d’una variable quantitativa contínua entre tres o més grups independents.

És exactament un cas particular del model lineal general (\(lm\)), on la variable predictora és un factor categòric:

\[Y = \beta_0 + \beta_1 X_{\text{grup 2}} + \beta_2 X_{\text{grup 3}} + \dots + \varepsilon\]

Analitzarem si existien diferències significatives en l’edat mitjana segons la classe de passatge (Primera, Segona o Tercera classe) a la cohort del Titanic.

1.4.1 Preparació de les dades

Codi 
titanic <- read_csv("titanic.csv")

dades_titanic <- titanic %>% 
  select(PassengerId, Survived, Pclass, Sex, Age) %>% 
  filter(!is.na(Age)) %>% 
  mutate(
    Pclass = factor(Pclass, levels = c(1, 2, 3), labels = c("Primera", "Segona", "Tercera")),
    Sex = factor(Sex, levels = c("female", "male"), labels = c("Dona", "Home")),
    Survived = factor(Survived, levels = c(0, 1), labels = c("Mort", "Viu"))
  )

1.4.2 Exploració visual prèvia

Codi 
ggplot(dades_titanic, aes(x = Pclass, y = Age, fill = Pclass)) +
  geom_boxplot(alpha = 0.7, outlier.alpha = 0.4) +
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3.5, color = "darkred") +
  scale_fill_brewer(palette = "Blues") +
  labs(title = "Distribució de l'edat segons la classe de passatge",
       subtitle = "Els rombes vermells indiquen la mitjana de cada grup",
       x = "Classe de passatge",
       y = "Edat (anys)") +
  theme_classic() +
  theme(legend.position = "none")
Figura 12

1.4.3 Ajust del model lineal i taula ANOVA

Codi 
model_anova <- lm(Age ~ Pclass, data = dades_titanic)
anova(model_anova)

1.4.3.1 Interpretació del contrast F:

  • L’estadístic \(F\) avalua la hipòtesi nul·la (\(H_0\)): totes les mitjanes poblacionals són iguals (\(\mu_{\text{Primera}} = \mu_{\text{Segona}} = \mu_{\text{Tercera}}\)).
  • L’obtenció d’un valor \(F = 57.44\) (\(p < 2.2 \times 10^{-16}\)) permet rebutjar la hipòtesi nul·la: existeixen diferències altament significatives en l’edat mitjana entre les classes de bitllet.

1.4.4 Verificació de supòsits del model ANOVA

Supòsit Què implica Com es comprova
Linealitat Garantida directament en ser un predictor categòric —
Independència Observacions no aparellades Disseny de l’estudi
Normalitat dels residus Els errors al voltant de les mitjanes segueixen una distribució normal Gràfic Q-Q i shapiro.test()
Homocedasticitat Variància de l’edat homogènia entre els tres grups Residus vs. ajustats i car::leveneTest()

1.4.4.1 Comprovació de la normalitat dels residus

Codi 
shapiro.test(residuals(model_anova))

    Shapiro-Wilk normality test

data:  residuals(model_anova)
W = 0.99146, p-value = 0.0003856
Codi 
plot(model_anova, which = 2)
Figura 13

Criteri clínic davant mostres grans (\(n > 50\)):
El test de Shapiro-Wilk resulta estadísticament significatiu (\(p < 0.001\)) a causa de l’elevada potència de la mostra (\(n = 714\)). Com que al gràfic Q-Q la immensa majoria d’observacions s’alinea fidelment sobre la diagonal teòrica, el model és prou robust i es pot continuar amb la interpretació paramètrica indicant aquesta lleu desviació.

1.4.4.2 Comprovació de l’homocedasticitat

Codi 
library(car)
leveneTest(Age ~ Pclass, data = dades_titanic)

1.4.5 Alternativa no paramètrica: Test de Kruskal-Wallis

Codi 
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-16

1.4.6 Interpretació detallada dels coeficients del model (summary)

Codi 
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-16
  • Intercept (\(\beta_0 = 38.23\) anys): Correspon exactament a l’edat mitjana dels passatgers de Primera classe.
  • PclassSegona (\(\beta_1 = -8.36\) anys; \(p < 0.001\)): Els passatgers de Segona classe tenen, de mitjana, 8.36 anys menys que els de Primera classe.
  • PclassTercera (\(\beta_2 = -13.09\) anys; \(p < 0.001\)): Els passatgers de Tercera classe tenen, de mitjana, 13.09 anys menys que els de Primera classe.

1.4.7 Comparacions múltiples post-hoc

1.4.7.1 1. Mètode paramètric: Test HSD de Tukey

Codi 
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.0004884

1.4.7.2 2. Mètode no paramètric: Test de Dunn

Codi 
library(FSA)
dunnTest(Age ~ Pclass, data = dades_titanic, method = "holm")
         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-05

1.5 Regressió logística binària

Quan la variable dependent clínica és binària o dicotòmica (viu/mort, recidiva sí/no, èxit terapèutic sí/no), la regressió lineal clàssica no és aplicable.

La regressió logística transforma la probabilitat d’un esdeveniment mitjançant la funció logit, assegurant que les prediccions es mantinguin sempre estrictament dins l’interval \([0, 1]\).

Model Variable resposta (\(Y\)) Resultat estimat (\(\hat{Y}\)) Aplicació clínica típica
Regressió lineal Quantitativa contínua Valor mitjà Pressió arterial sistòlica, colesterol LDL
Regressió logística Binària (0/1) Probabilitat (0 a 1) Supervivència a 30 dies, risc d’ingrés a UCI

1.5.1 La funció logística i l’escala Logit

\[\text{logit}(p) = \ln\left(\frac{p}{1 - p}\right) = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \dots + \beta_k X_k\]

Codi 
curva_logistica <- tibble(
  x = seq(-6, 6, length.out = 200),
  p = 1 / (1 + exp(-x))
)

ggplot(curva_logistica, aes(x = x, y = p)) +
  geom_line(color = "#1976D2", linewidth = 1.5) +
  geom_hline(yintercept = c(0, 1), linetype = "dotted", color = "gray50") +
  labs(title = "Funció sigmoide de la regressió logística",
       subtitle = "La probabilitat queda confinada asimptòticament entre 0 i 1",
       x = "Combinació lineal de predictors (logit)",
       y = "Probabilitat de l'esdeveniment (p)") +
  theme_classic()
Figura 14

1.5.2 Ajust del model logístic amb glm()

Codi 
model_log <- glm(Survived ~ Sex + Pclass + Age, 
                 data = dades_titanic, 
                 family = binomial(link = "logit"))

summary(model_log)

Call:
glm(formula = Survived ~ Sex + Pclass + Age, family = binomial(link = "logit"), 
    data = dades_titanic)

Coefficients:
               Estimate Std. Error z value Pr(>|z|)    
(Intercept)    3.777013   0.401123   9.416  < 2e-16 ***
SexHome       -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: 5

1.5.3 De l’escala Log-Odds a les Odds Ratios (OR)

\[\text{OR} = e^{\beta}\]

  • \(\text{OR} = 1\): El factor no modifica la probabilitat de l’esdeveniment.
  • \(\text{OR} > 1\): El factor s’associa a un major risc o probabilitat de l’esdeveniment.
  • \(\text{OR} < 1\): El factor actua com a protector, disminuint la probabilitat.
Codi 
round(exp(coef(model_log)), 3)
  (Intercept)       SexHome  PclassSegona PclassTercera           Age 
       43.685         0.080         0.270         0.076         0.964
Codi 
round(exp(confint(model_log)), 3)
               2.5 % 97.5 %
(Intercept)   20.379 98.386
SexHome        0.053  0.119
PclassSegona   0.155  0.462
PclassTercera  0.043  0.130
Age            0.949  0.978

1.5.3.1 Interpretació clínica dels coeficients:

  • SexHome (\(\text{OR} = 0.080\), IC 95%: \(0.053 - 0.119\); \(p < 0.001\)): Ajustant per classe i edat, els homes presenten una probabilitat de sobreviure un 92% inferior a la de les dones (\(\text{OR} = 0.08\)).
  • PclassSegona (\(\text{OR} = 0.270\), IC 95%: \(0.155 - 0.462\); \(p < 0.001\)): Viatjar a Segona classe s’associa a una reducció del 73% en els odds de supervivència respecte a la Primera classe de referència.
  • PclassTercera (\(\text{OR} = 0.076\), IC 95%: \(0.043 - 0.130\); \(p < 0.001\)): Viatjar a Tercera classe redueix els odds de sobreviure un 92.4% en comparació amb Primera classe.
  • Age (\(\text{OR} = 0.964\), IC 95%: \(0.949 - 0.978\); \(p < 0.001\)): Per cada any addicional d’edat, els odds de supervivència disminueixen aproximadament un 3.6% (\(\text{OR} = 0.964\)).

1.5.4 Verificació de supòsits del model logístic

  1. Independència de les observacions: No hi ha mesures repetides sobre el mateix pacient sense models mixts.
  2. Linealitat del logit: Els predictors continus han de guardar una relació lineal amb el logit de la resposta.
  3. Absència de multicol·linealitat: Els predictors no han d’estar excessivament correlacionats entre ells.

1.5.4.1 Linealitat del logit per a l’edat

Codi 
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.25, color = "gray30") +
  geom_smooth(method = "loess", se = FALSE, color = "#D32F2F", linewidth = 1.2) +
  labs(title = "Avaluació de la linealitat del logit respecte a l'edat",
       x = "Edat (anys)", 
       y = "Logit predit") +
  theme_classic()
Figura 15

1.5.4.2 Absència de col·linealitat (VIF)

Codi 
library(car)
vif(model_log)
           GVIF Df GVIF^(1/(2*Df))
Sex    1.072870  1        1.035794
Pclass 1.416574  2        1.090962
Age    1.333112  1        1.154605

1.5.5 Capacitat discriminativa: Corba ROC i Àrea Sota la Corba (AUC)

Codi 
library(pROC)

roc_obj <- roc(dades_titanic$Survived, dades_titanic$prob_pred, levels = c("Mort", "Viu"))

plot(roc_obj, 
     col = "#1976D2", 
     lwd = 2.5, 
     main = "Corba ROC: Model predictiu de supervivència",
     xlab = "1 - Especificitat", 
     ylab = "Sensibilitat")
abline(a = 0, b = 1, lty = 2, col = "gray60")

auc(roc_obj)
Area under the curve: 0.8523
Codi 
ci.auc(roc_obj)
95% CI: 0.8219-0.8827 (DeLong)
Figura 16

1.5.5.1 Criteri clínic d’interpretació de l’AUC:

  • \(\text{AUC} = 0.50\): Capacitat nul·la (equivalent a l’atzar).
  • \(\text{AUC } 0.70 - 0.79\): Discriminació acceptable.
  • \(\text{AUC } 0.80 - 0.89\): Molt bona discriminació (el nostre model obté 0.852).
  • \(\text{AUC } \ge 0.90\): Discriminació excel·lent (vigilar possible sobreajust).

1.5.6 Taules de regressió per a publicació amb gtsummary

Codi 
model_log %>% 
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      Sex    ~ "Sexe",
      Pclass ~ "Classe de passatge",
      Age    ~ "Edat (anys)"
    )
  ) %>% 
  bold_p(t = 0.05) %>% 
  bold_labels()
Characteristic OR 95% CI p-value
Sexe


    Dona — —
    Home 0.08 0.05, 0.12 <0.001
Classe de passatge


    Primera — —
    Segona 0.27 0.16, 0.46 <0.001
    Tercera 0.08 0.04, 0.13 <0.001
Edat (anys) 0.96 0.95, 0.98 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

1.6 Síntesi: Comparativa dels models estudiats

Model Funció a R Resposta (\(Y\)) Coeficients clau Supòsits essencials
Regressió lineal simple lm(y ~ x) Quantitativa contínua Pendent (\(\beta_1\), canvi en unitats de \(Y\)) Linealitat, normalitat residus, homocedasticitat
Regressió lineal múltiple lm(y ~ x1 + x2) Quantitativa contínua Efecte ajustat per confusors Supòsits lineals + absència de col·linealitat (VIF)
ANOVA d’un factor lm(y ~ factor) Quantitativa contínua Diferència de mitjanes vs. grup de referència Normalitat dels residus, homocedasticitat
Regressió logística glm(..., family = binomial) Categòrica binària (0/1) Odds Ratios (\(\text{OR} = e^\beta\)) Linealitat del logit, independència, absència de col·linealitat

1.7 Exercicis pràctics del Capítol 4

  1. Regressió lineal simple: Utilitzant el conjunt de dades iris, ajusta un model lineal per predir la longitud del pètal (Petal.Length) a partir de la seva amplada (Petal.Width). Interpreta el pendent obtingut, el contrast d’hipòtesi i el coeficient \(R^2\).
  2. Diagnòstic de supòsits: Genera la bateria de 4 gràfics diagnòstics del model anterior (plot(model)). Avalua visualment la linealitat, l’homocedasticitat i comprova formalment la normalitat dels residus amb shapiro.test().
  3. Regressió lineal múltiple: Ajusta un model múltiple per estimar Petal.Length a partir de Petal.Width i Sepal.Length. Avalua si tots dos predictors aporten informació independent i comprova que no hi hagi col·linealitat mitjançant car::vif().
  4. ANOVA com a model lineal: Comprova si existeixen diferències significatives en la longitud del pètal (Petal.Length) entre les tres espècies (Species) d’iris mitjançant lm() i anova(). Realitza les comparacions post-hoc amb TukeyHSD().
  5. Regressió logística binària: A la cohort del Titanic (data/titanic.csv), crea una variable binària per identificar si el bitllet va costar més de 30 dòlars (tarifa_alta = ifelse(Fare > 30, 1, 0)). Ajusta un model logístic amb glm() per predir tarifa_alta en funció de la classe (Pclass) i el sexe (Sex). Obtén i interpreta les Odds Ratios (OR) amb els seus intervals de confiança del 95%.
  6. Capacitat discriminativa (ROC/AUC): Extreu les probabilitats predites del model logístic anterior, dibuixa la corba ROC amb el paquet pROC i calcula l’Àrea Sota la Corba (AUC) amb el seu IC 95%.
  7. Taula de publicació: Genera una taula elegant i a punt per a publicació amb gtsummary::tbl_regression() del model logístic anterior, mostrant les OR exponenciades i els noms de les variables personalitzats.

1.7.1 Solucions als exercicis (fes clic per desplegar)

Codi 
library(tidyverse)
library(car)
library(pROC)
library(gtsummary)

# Exercici 1: Regressió lineal simple amb iris
model_iris_simple <- lm(Petal.Length ~ Petal.Width, data = iris)
summary(model_iris_simple)

Call:
lm(formula = Petal.Length ~ Petal.Width, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.33542 -0.30347 -0.02955  0.25776  1.39453 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.08356    0.07297   14.85   <2e-16 ***
Petal.Width  2.22994    0.05140   43.39   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4782 on 148 degrees of freedom
Multiple R-squared:  0.9271,    Adjusted R-squared:  0.9266 
F-statistic:  1882 on 1 and 148 DF,  p-value: < 2.2e-16
Codi 
# Exercici 2: Diagnòstic de supòsits
par(mfrow = c(2, 2))
plot(model_iris_simple)

Codi 
par(mfrow = c(1, 1))

shapiro.test(residuals(model_iris_simple))

    Shapiro-Wilk normality test

data:  residuals(model_iris_simple)
W = 0.99011, p-value = 0.3753
Codi 
# Exercici 3: Regressió lineal múltiple i VIF
model_iris_mult <- lm(Petal.Length ~ Petal.Width + Sepal.Length, data = iris)
summary(model_iris_mult)

Call:
lm(formula = Petal.Length ~ Petal.Width + Sepal.Length, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.15506 -0.21920 -0.02115  0.25986  1.35204 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.50714    0.33696  -4.473 1.54e-05 ***
Petal.Width   1.74810    0.07533  23.205  < 2e-16 ***
Sepal.Length  0.54226    0.06934   7.820 9.41e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4032 on 147 degrees of freedom
Multiple R-squared:  0.9485,    Adjusted R-squared:  0.9478 
F-statistic:  1354 on 2 and 147 DF,  p-value: < 2.2e-16
Codi 
vif(model_iris_mult)
 Petal.Width Sepal.Length 
    3.021401     3.021401
Codi 
# Exercici 4: ANOVA com a model lineal i Post-hoc
model_anova_iris <- lm(Petal.Length ~ Species, data = iris)
anova(model_anova_iris)
Codi 
TukeyHSD(aov(model_anova_iris))
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = model_anova_iris)

$Species
                      diff     lwr     upr p adj
versicolor-setosa    2.798 2.59422 3.00178     0
virginica-setosa     4.090 3.88622 4.29378     0
virginica-versicolor 1.292 1.08822 1.49578     0
Codi 
# Exercici 5: Regressió logística amb el Titanic
titanic <- read_csv("titanic.csv")

dades_ex5 <- titanic %>% 
  filter(!is.na(Fare), !is.na(Pclass), !is.na(Sex)) %>% 
  mutate(
    tarifa_alta = ifelse(Fare > 30, 1, 0),
    Pclass = factor(Pclass, levels = c(1, 2, 3), labels = c("Primera", "Segona", "Tercera")),
    Sex = factor(Sex, levels = c("male", "female"), labels = c("Home", "Dona"))
  )

model_log_ex5 <- glm(tarifa_alta ~ Pclass + Sex, data = dades_ex5, family = binomial)

round(exp(coef(model_log_ex5)), 3)
  (Intercept)  PclassSegona PclassTercera       SexDona 
        2.667         0.038         0.026         2.192
Codi 
round(exp(confint(model_log_ex5)), 3)
              2.5 % 97.5 %
(Intercept)   1.886  3.836
PclassSegona  0.022  0.066
PclassTercera 0.016  0.040
SexDona       1.448  3.338
Codi 
# Exercici 6: Corba ROC i AUC
dades_ex5$probabilitat <- predict(model_log_ex5, type = "response")

roc_ex6 <- roc(dades_ex5$tarifa_alta, dades_ex5$probabilitat)
plot(roc_ex6, col = "#1976D2", lwd = 2, main = "Corba ROC - Predicció de Tarifa > 30$")

Codi 
auc(roc_ex6)
Area under the curve: 0.8531
Codi 
ci.auc(roc_ex6)
95% CI: 0.8209-0.8853 (DeLong)
Codi 
# Exercici 7: Taula de publicació amb gtsummary
model_log_ex5 %>% 
  tbl_regression(
    exponentiate = TRUE,
    label = list(
      Pclass ~ "Classe de bitllet",
      Sex    ~ "Sexe"
    )
  ) %>% 
  bold_p(t = 0.05) %>% 
  bold_labels()
Characteristic OR 95% CI p-value
Classe de bitllet


    Primera — —
    Segona 0.04 0.02, 0.07 <0.001
    Tercera 0.03 0.02, 0.04 <0.001
Sexe


    Home — —
    Dona 2.19 1.45, 3.34 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio