Bioestadística Fundamental y Estadística Fundamental para las Ciencias de la Salud
Authors
Affiliations
Jose Miguel Leon Puentes
Universidad Nacional de Colombia
Johann David Sanchez Niño
Andrés Felipe Rache
Laura Camila Matheus Barrios
Brayan Stiven Martínez Rodríguez
Departamento de Estadística
Ejercicio 1. Edad y presión arterial sistólica
Solución
# a) Datos y ajuste del modeloedad <-c(25, 31, 38, 42, 29, 54, 60, 47, 33, 58, 45, 22, 36, 51, 63)pas <-c(125, 121, 135, 137, 122, 134, 151, 137, 130, 150, 138, 115, 131, 145, 143)modelo1 <-lm(pas ~ edad)summary(modelo1)
Call:
lm(formula = pas ~ edad)
Residuals:
Min 1Q Median 3Q Max
-8.894 -3.439 1.724 3.562 4.312
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 103.19011 4.09924 25.173 2.05e-12 ***
edad 0.73525 0.09293 7.912 2.52e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.544 on 13 degrees of freedom
Multiple R-squared: 0.828, Adjusted R-squared: 0.8148
F-statistic: 62.6 on 1 and 13 DF, p-value: 2.523e-06
edad y pas son los dos vectores numéricos con la información de los 15 pacientes. La función lm(pas ~ edad) ajusta un modelo lineal (“linear model”) en el que, a la izquierda del símbolo ~, se indica la variable respuesta (pas) y, a la derecha, la variable explicativa (edad); R estima los coeficientes \(\hat\beta_0\) (intercepto) y \(\hat\beta_1\) (pendiente) por mínimos cuadrados ordinarios. summary() muestra un resumen completo del modelo: los coeficientes estimados (Estimate), su error estándar (Std. Error), el estadístico \(t = \text{Estimate}/\text{Std. Error}\) y su \(p\)-valor asociado (columna Pr(>|t|)), además del \(R^2\), el \(R^2\) ajustado y el estadístico F global del modelo.
En este caso, \(\hat\beta_0 = 103.19\) y \(\hat\beta_1 = 0.7353\), con \(R^2 = 0.828\).
b) Interpretación de los coeficientes
El intercepto (\(\hat\beta_0 = 103.19\) mmHg) representa la presión arterial sistólica esperada cuando la edad es igual a cero; en este contexto no tiene una interpretación clínica directa (nadie tiene edad cero en la muestra), por lo que se interpreta principalmente como un valor de referencia matemático del modelo. La pendiente (\(\hat\beta_1 = 0.7353\)) indica que, por cada año adicional de edad, la presión arterial sistólica esperada aumenta en promedio 0.7353 mmHg, manteniendo el resto de condiciones constantes.
c) Significancia de la pendiente
# c) p-valor asociado a la pendientecoef(summary(modelo1))["edad", ]
Estimate Std. Error t value Pr(>|t|)
7.352498e-01 9.292798e-02 7.912040e+00 2.523376e-06
coef(summary(modelo1)) extrae la tabla de coeficientes del resumen del modelo en forma de matriz; al indexarla con ["edad", ] obtenemos únicamente la fila correspondiente a la pendiente de edad, con su estimación, error estándar, estadístico t y p-valor. El p-valor obtenido es \(2.52\times10^{-6}\), muy inferior a \(\alpha = 0.05\). Por lo tanto, se rechaza \(H_0: \beta_1 = 0\) y se concluye que existe evidencia estadísticamente significativa de que la edad tiene un efecto lineal sobre la presión arterial sistólica.
d) Coeficiente de determinación
# d) R^2summary(modelo1)$r.squared
[1] 0.8280432
summary(modelo1)$r.squared accede directamente al valor de \(R^2\) almacenado dentro del objeto que produce summary(). Se obtiene \(R^2 = 0.828\), lo que significa que el 82.8 % de la variabilidad observada en la presión arterial sistólica es explicada linealmente por la edad; el 17.2 % restante se debe a otros factores no incluidos en el modelo o a variabilidad aleatoria.
# Predicción con predict()predict(modelo1, newdata =data.frame(edad =55))
1
143.6288
coef(modelo1) devuelve el vector de coeficientes estimados; coef(modelo1)[1] es el intercepto y coef(modelo1)[2] es la pendiente. Calculamos manualmente \(\hat{y} = \hat\beta_0 + \hat\beta_1 \times 55\). Alternativamente, predict(modelo1, newdata = data.frame(edad = 55)) le pide al modelo que evalúe la ecuación ajustada en un nuevo valor de edad; el argumento newdata debe ser un data frame cuyas columnas tengan el mismo nombre que las variables explicativas usadas en lm(). Ambos caminos coinciden: 143.63 mmHg es la presión arterial sistólica esperada para una persona de 55 años, según el modelo ajustado.
Ejercicio 2. Dosis del fármaco y reducción del colesterol
Solución
# a) Datos y ajuste del modelodosis <-c(5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60)reduccion <-c(10.5, 10.7, 12.2, 11.2, 20.3, 19.3, 23.0, 21.2, 35.8, 33.2, 43.7, 46.1)modelo2 <-lm(reduccion ~ dosis)summary(modelo2)
Call:
lm(formula = reduccion ~ dosis)
Residuals:
Min 1Q Median 3Q Max
-7.795 -2.702 0.752 3.476 5.124
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.00152 2.60500 0.768 0.46
dosis 0.67483 0.07079 9.533 2.46e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4.233 on 10 degrees of freedom
Multiple R-squared: 0.9009, Adjusted R-squared: 0.891
F-statistic: 90.87 on 1 and 10 DF, p-value: 2.46e-06
Igual que en el ejercicio anterior, lm(reduccion ~ dosis) ajusta el modelo de regresión lineal simple con reduccion como respuesta y dosis como predictor.
b) Significancia del intercepto
El resumen muestra \(\hat\beta_0 = 2.0015\) con un \(p\)-valor de \(0.46\) para el intercepto. Como \(0.46 > 0.05\), no hay evidencia suficiente para rechazar \(H_0: \beta_0 = 0\); es decir, el intercepto no es significativamente distinto de cero. Esto tiene sentido clínico: con una dosis de 0 mg (ausencia de tratamiento) es razonable esperar una reducción de colesterol cercana a 0 mg/dL.
c) Interpretación de la pendiente
La pendiente estimada es \(\hat\beta_1 = 0.6748\) (\(p = 2.46\times10^{-6}\), altamente significativa). Esto significa que, por cada miligramo adicional de dosis del fármaco, la reducción esperada de colesterol aumenta en promedio 0.6748 mg/dL.
d) Análisis de residuales
# d) Residuales vs valores ajustadosplot(fitted(modelo2), resid(modelo2),xlab ="Valores ajustados", ylab ="Residuales",main ="Residuales vs. valores ajustados - Modelo 2", pch =19, col ="#457B9D")abline(h =0, col ="red", lty =2)
fitted(modelo2) devuelve los valores predichos (\(\hat{y}_i\)) por el modelo para cada observación de la muestra, y resid(modelo2) devuelve los residuales (\(e_i = y_i - \hat{y}_i\)). Al graficar los residuales contra los valores ajustados esperamos, bajo el supuesto de homocedasticidad (varianza constante de los errores), una nube de puntos dispersa aleatoriamente alrededor de la línea horizontal h = 0 (agregada con abline(h = 0, ...)), sin ningún patrón sistemático. En este caso, la dispersión de los residuales tiende a aumentar levemente a medida que crecen los valores ajustados (los residuales son más pequeños para dosis bajas y más grandes en valor absoluto para dosis altas), lo cual es un indicio visual leve de heterocedasticidad. En la práctica, ante esta señal se recomendaría explorar una transformación de la variable respuesta (por ejemplo, \(\sqrt{\text{reduccion}}\) o \(\log(\text{reduccion})\)) o usar errores estándar robustos a heterocedasticidad.
Ejercicio 3. Capacidad pulmonar en función de estatura y edad
Call:
lm(formula = fev ~ altura * edad3)
Residuals:
Min 1Q Median 3Q Max
-0.45834 -0.13156 -0.01726 0.13138 0.78378
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.642499 5.423925 -0.118 0.907
altura 0.023031 0.032917 0.700 0.494
edad3 -0.075276 0.383532 -0.196 0.847
altura:edad3 0.001529 0.002337 0.654 0.522
Residual standard error: 0.2882 on 16 degrees of freedom
Multiple R-squared: 0.8236, Adjusted R-squared: 0.7905
F-statistic: 24.9 on 3 and 16 DF, p-value: 2.873e-06
La notación altura * edad3 dentro de lm() es una forma abreviada de escribir altura + edad3 + altura:edad3: incluye los dos efectos principales (altura, edad3) y su término de interacción (altura:edad3). El coeficiente de interacción mide si el efecto de la estatura sobre el FEV cambia según la edad (o, equivalentemente, si el efecto de la edad cambia según la estatura).
b) Significancia de la interacción
En la tabla de summary(), el coeficiente altura:edad3 es \(0.001529\) con \(p = 0.522\), muy por encima de \(\alpha = 0.05\): individualmente no es significativo. Confirmamos con una prueba F comparando el modelo reducido (sin interacción) contra el modelo completo:
# Comparación de modelos con anova()modelo3_red <-lm(fev ~ altura + edad3)anova(modelo3_red, modelo3_int)
Analysis of Variance Table
Model 1: fev ~ altura + edad3
Model 2: fev ~ altura * edad3
Res.Df RSS Df Sum of Sq F Pr(>F)
1 17 1.3649
2 16 1.3294 1 0.035534 0.4277 0.5224
anova(modelo3_red, modelo3_int) compara dos modelos anidados (el primero es un caso particular del segundo, fijando el coeficiente de interacción en cero) mediante una prueba F que evalúa si la reducción en la suma de cuadrados residual (Sum Sq) al pasar del modelo reducido al completo es estadísticamente significativa. Se obtiene \(F = 0.4277\) con \(p = 0.5224\). Como \(p > 0.05\), no se rechaza \(H_0\): la interacción no aporta una mejora significativa al modelo.
c) Modelo final
Dado que la interacción no es significativa, se opta por el modelo más simple (sin interacción):
# c) Modelo final sin interacciónsummary(modelo3_red)
Call:
lm(formula = fev ~ altura + edad3)
Residuals:
Min 1Q Median 3Q Max
-0.45968 -0.11908 -0.01659 0.12374 0.77665
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.100897 1.184621 -3.462 0.00298 **
altura 0.044107 0.006583 6.700 3.73e-06 ***
edad3 0.175015 0.024399 7.173 1.56e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2834 on 17 degrees of freedom
Multiple R-squared: 0.8189, Adjusted R-squared: 0.7976
F-statistic: 38.43 on 2 and 17 DF, p-value: 4.925e-07
Con este modelo, \(\hat\beta_{altura} = 0.04411\) (\(p = 3.73\times10^{-6}\)) y \(\hat\beta_{edad3} = 0.17502\) (\(p = 1.56\times10^{-6}\)), ambos altamente significativos. Interpretación (excepto el intercepto): manteniendo la edad constante, por cada centímetro adicional de estatura, el FEV esperado aumenta en promedio 0.0441 litros; manteniendo la estatura constante, por cada año adicional de edad, el FEV esperado aumenta en promedio 0.175 litros.
d) Predicción
# d) Predicción para altura = 165 cm, edad = 13 añospredict(modelo3_red, newdata =data.frame(altura =165, edad3 =13))
1
5.451941
Usando el modelo final (sin interacción), la capacidad pulmonar esperada para un adolescente de 165 cm y 13 años es de 5.45 litros aproximadamente.
Ejercicio 4. Pérdida de peso según horas de ejercicio y tipo de dieta
Solución
# a) Construcción del data framehoras <-rep(1:8, times =2)dieta <-factor(c(rep("A", 8), rep("B", 8)))perdida <-c(1.18, 1.16, 2.01, 1.64, 3.10, 2.44, 2.68, 3.68,2.04, 2.40, 2.80, 4.27, 5.47, 4.42, 6.49, 6.41)datos4 <-data.frame(horas, dieta, perdida)head(datos4)
horas dieta perdida
1 1 A 1.18
2 2 A 1.16
3 3 A 2.01
4 4 A 1.64
5 5 A 3.10
6 6 A 2.44
rep(1:8, times = 2) repite la secuencia de horas (1 a 8) dos veces, una por cada dieta. factor(c(rep("A", 8), rep("B", 8))) crea la variable categórica dieta con 8 repeticiones de "A" seguidas de 8 de "B", y factor() le indica a R que esta variable de texto debe tratarse como una variable cualitativa (con niveles), no como texto libre. data.frame() combina los tres vectores en una tabla de datos; head() muestra las primeras 6 filas para verificar que la construcción fue correcta.
b) Modelo con interacción
# b) Ajuste del modelo con interacción horas:dietamodelo4 <-lm(perdida ~ horas * dieta, data = datos4)summary(modelo4)
Call:
lm(formula = perdida ~ horas * dieta, data = datos4)
Residuals:
Min 1Q Median 3Q Max
-0.88714 -0.31646 -0.04196 0.29310 0.84262
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.74429 0.41012 1.815 0.09461 .
horas 0.33155 0.08122 4.082 0.00152 **
dietaB 0.48429 0.57999 0.835 0.42004
horas:dietaB 0.34821 0.11486 3.032 0.01043 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5263 on 12 degrees of freedom
Multiple R-squared: 0.9248, Adjusted R-squared: 0.9059
F-statistic: 49.16 on 3 and 12 DF, p-value: 5.149e-07
Al usar data = datos4, lm() toma las variables perdida, horas y dieta directamente de las columnas del data frame. Como dieta es un factor con niveles "A" y "B", R toma automáticamente "A" como categoría de referencia (por ser la primera en orden alfabético) y crea una variable indicadora dietaB (que vale 1 si la dieta es B, y 0 si es A). El término horas:dietaB es la interacción entre las horas de ejercicio y el hecho de estar en la dieta B.
El coeficiente de interacción es \(\hat\beta_{horas:dietaB} = 0.3482\) con \(p = 0.0104 < 0.05\): es significativo. Esto indica que el efecto de las horas de ejercicio sobre la pérdida de peso sí depende del tipo de dieta: en la dieta B, cada hora adicional de ejercicio se asocia con 0.3482 kg adicionales de pérdida de peso, por encima del efecto que ya tienen las horas de ejercicio en la dieta A.
c) Interpretación de la interacción
Con dieta A (categoría de referencia), la pendiente de horas es directamente \(\hat\beta_{horas} = 0.3316\) kg por hora de ejercicio. Con dieta B, la pendiente efectiva de horas es \(\hat\beta_{horas} + \hat\beta_{horas:dietaB} = 0.3316 + 0.3482 = 0.6797\) kg por hora de ejercicio, es decir, más del doble que con la dieta A. En otras palabras, cada hora de ejercicio “rinde” más kilogramos perdidos bajo la dieta B que bajo la dieta A.
d) Interpretación de dietaB
El coeficiente \(\hat\beta_{dietaB} = 0.4843\) (\(p = 0.42\), no significativo) representa la diferencia esperada en pérdida de peso entre la dieta B y la dieta A cuando horas = 0; al no ser significativo (y al no tener sentido evaluarlo en horas = 0, un valor fuera del rango observado y sin relevancia práctica), no se interpreta como una diferencia relevante entre dietas “en reposo”. La comparación relevante entre dietas es, más bien, la de las pendientes vista en el punto c).
e) Predicción
# e) Predicción con predict()predict(modelo4, newdata =data.frame(horas =6, dieta =factor("B", levels =c("A", "B"))))
1
5.307143
Se construye un data frame de una sola fila con horas = 6 y dieta = "B" (usando factor(..., levels = c("A","B")) para asegurar que el nuevo dato use los mismos niveles del factor que el modelo original). La pérdida de peso esperada para una persona que hace 6 horas de ejercicio semanal con la dieta B es de 5.31 kg aproximadamente.
Ejercicio 5. Colesterol en función de edad, IMC y ejercicio
Solución
# a) Datos (n = 25, 4 variables), cargados desde el archivo externo# "mod7_E05_datos.csv" (columnas: "edad", "imc", "ejercicio_h", "colesterol").datos5 <-read.csv("mod7_E05_datos.csv")edad5 <- datos5$edadimc5 <- datos5$imcejercicio5 <- datos5$ejercicio_hcolesterol <- datos5$colesterolmodelo5 <-lm(colesterol ~ edad5 + imc5 + ejercicio5)summary(modelo5)
Call:
lm(formula = colesterol ~ edad5 + imc5 + ejercicio5)
Residuals:
Min 1Q Median 3Q Max
-20.1951 -6.9816 0.9971 4.8009 20.3284
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 161.6791 14.6493 11.037 3.35e-10 ***
edad5 0.9987 0.2308 4.327 0.000297 ***
imc5 2.2155 0.6215 3.565 0.001831 **
ejercicio5 -4.5553 1.2380 -3.680 0.001394 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 9.316 on 21 degrees of freedom
Multiple R-squared: 0.828, Adjusted R-squared: 0.8034
F-statistic: 33.69 on 3 and 21 DF, p-value: 3.274e-08
read.csv("mod7_E05_datos.csv") lee el archivo CSV con las 25 observaciones y sus 4 variables (debe estar guardado en la misma carpeta que este documento) y lo carga como un data frame; cada $columna extrae una variable como vector numérico, equivalente a escribir los 25 valores de cada una manualmente con c(...).
El modelo incluye tres predictores sin interacciones. Los tres coeficientes resultan significativos al 5 %: \(\hat\beta_{edad5} = 0.999\) (\(p < 0.001\)), \(\hat\beta_{imc5} = 2.216\) (\(p = 0.0018\)) y \(\hat\beta_{ejercicio5} = -4.555\) (\(p = 0.0014\)). Interpretación (manteniendo las demás variables constantes): por cada año adicional de edad, el colesterol esperado aumenta 0.999 mg/dL; por cada unidad adicional de IMC, aumenta 2.216 mg/dL; y por cada hora adicional de ejercicio semanal, el colesterol esperado disminuye 4.555 mg/dL. El modelo explica \(R^2 = 0.828\) de la variabilidad del colesterol.
b) Colinealidad: correlación y VIF
# b) Correlación entre las variables explicativascor(edad5, imc5)
[1] 0.5145222
cor(edad5, imc5) calcula el coeficiente de correlación de Pearson entre las dos variables. Se obtiene \(r = 0.51\), una correlación positiva moderada: las personas de mayor edad tienden a tener, en promedio, un IMC algo más alto en esta muestra. Una correlación moderada como esta no necesariamente es un problema, pero conviene cuantificar su impacto sobre el modelo con el VIF.
# VIF calculado "a mano" (sin paquetes adicionales)# El VIF de una variable X_j se obtiene ajustando X_j en función de las demás# variables explicativas y usando el R^2 de esa regresión auxiliar:# VIF_j = 1 / (1 - R^2_j)vif_manual <-function(modelo) { X <-model.matrix(modelo)[, -1, drop =FALSE] # quita la columna del interceptosapply(colnames(X), function(v) { otras <-setdiff(colnames(X), v) aux <-lm(X[, v] ~ X[, otras])1/ (1-summary(aux)$r.squared) })}vif_manual(modelo5)
edad5 imc5 ejercicio5
1.394734 1.399176 1.111717
Esta función replica manualmente lo que hace el Factor de Inflación de la Varianza: para cada predictor \(X_j\), se ajusta una regresión auxiliar de \(X_j\) contra el resto de predictores del modelo, se obtiene su \(R^2_j\), y se calcula \(VIF_j = 1/(1-R^2_j)\). Cuanto mayor es la capacidad de las otras variables para “explicar” a \(X_j\) (mayor \(R^2_j\)), mayor es el VIF y más colineal es esa variable con las demás. model.matrix(modelo) construye la matriz de diseño del modelo (incluye el intercepto, que se descarta con [, -1]); setdiff() calcula el conjunto de columnas restantes. Si tuviera instalado el paquete car, el mismo resultado se obtendría de forma más directa con car::vif(modelo5).
Los VIF obtenidos son todos cercanos a 1 (aproximadamente 1.39, 1.40 y 1.11), muy por debajo del umbral de preocupación (5 o 10). Se concluye que, a pesar de la correlación moderada entre edad5 e imc5, no existe un problema de multicolinealidad relevante en este modelo: las tres variables aportan información suficientemente independiente entre sí.
c) Distancia de Cook
# c) Distancia de Cookcooks <-cooks.distance(modelo5)round(cooks, 4)
n <-length(colesterol)umbral_cook <-4/ numbral_cook
[1] 0.16
which(cooks > umbral_cook)
14 25
14 25
cooks.distance(modelo5) calcula, para cada observación, la distancia de Cook \(D_i\), una medida de cuánto cambiarían los coeficientes estimados del modelo si esa observación se eliminara de la muestra (combina el residual y el apalancamiento de cada punto). Un criterio práctico común es marcar como potencialmente influyente cualquier observación con \(D_i > 4/n\); aquí n <- length(colesterol) es el tamaño de muestra (25) y umbral_cook <- 4/n = 0.16. which(cooks > umbral_cook) devuelve los índices de las observaciones que superan ese umbral. Se obtienen las observaciones 14 (\(D_{14} \approx 0.207\)) y 25 (\(D_{25} \approx 0.235\)). Ninguna de las dos alcanza el umbral más estricto de \(D_i > 1\), así que no se trata de puntos que distorsionen gravemente el modelo, pero sí conviene revisarlos con más detalle.
d) Apalancamiento (leverage)
# d) Valores de apalancamiento (hat values)lev <-hatvalues(modelo5)round(lev, 4)
p <-3# numero de variables explicativas del modeloumbral_lev <-2* (p +1) / numbral_lev
[1] 0.32
which(lev > umbral_lev)
7 14
7 14
hatvalues(modelo5) calcula los elementos de la diagonal de la “matriz sombrero” (hat matrix), \(h_{ii}\), que miden qué tan alejada está cada observación del centro de la nube de valores de las variables explicativas (no dependen de la variable respuesta). El criterio habitual marca como de “alto apalancamiento” a las observaciones con \(h_{ii} > 2(p+1)/n\), donde \(p\) es el número de predictores (aquí \(p = 3\), por lo que el umbral es \(2(3+1)/25 = 0.32\)). Se identifican las observaciones 7 (\(h_{77} \approx 0.465\)) y 14 (\(h_{14,14} \approx 0.373\)) como puntos de alto apalancamiento: sus combinaciones de edad, IMC y horas de ejercicio son atípicas respecto al resto de la muestra.
e) Conclusión del diagnóstico
En conjunto, el modelo ajustado para el colesterol es razonablemente estable: no hay evidencia de multicolinealidad problemática entre los predictores (todos los VIF cercanos a 1), y aunque existen un par de observaciones con apalancamiento relativamente alto (7 y 14) y otro par con distancia de Cook algo mayor que el resto (14 y 25), ninguna alcanza los umbrales más estrictos que indicarían una observación claramente influyente y distorsionante (\(D_i > 1\)). La observación 14 aparece señalada tanto por apalancamiento alto como por distancia de Cook relativamente alta, por lo que sería la primera candidata a revisar (verificar que sus datos no correspondan a un error de digitación) antes de reportar el modelo como definitivo.
Referencias
Timbers, T., Campbell, T., & Lee, M. (2024). Data Science: A First Introduction. CRC Press. https://datasciencebook.ca/ (texto de referencia general del Módulo 7).
