# Activación de paquetes
library(haven)
library(foreign)
library(ggplot2)
library(moments)
library(knitr)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(olsrr)  # Métodos de selección por pasos
library(PerformanceAnalytics)  # Varias funciones.
# Lo utilizaremos para construir matriz de correlaciones, incluyendo gráficos
# de dispersión para análisis exploratorio de los datos

Taller evaluativo parte II.

Consideremos los datos del estudio “Heart and Estrogen/Progestin Study (HERS), un ensayo clínico de la terapia hormonal para la prevención de ataques cardíacos recurrentes y muertes entre 2.763 mujeres posmenopáusicas con enfermedad coronaria existente Grady et al. (1998).

Los investigadores desean conocer el efecto de diferentes factores sobre los niveles colesterol total (tchol1) al año de seguimiento. Los factores a analizar son: edad (age), raza (raceth:1=blanca, 2=Afroamericana, 3=otra raza), indice de masa corporal (BMI), puntaje subjetivo de autoreporte de salud (globrat: escala de 1 a 5, donde 5 indica un estado de salud “ideal”), tabaquismo (smoking), actividad física (physact: escala de 1 a 5, desde sedentario a actividad vigorosa), antecedente medico o enfermedad (medcond), diabetes (diabetes), colesterol basal (tchol), glucosa en sangre basal (glucose), uso de estatinas (statins), niveles basales de LDL y HDL y trigliceridos (TG), presión arterial sistólica (SBP) y diastólica (DBP).

# Cargar base de datos
library(haven)
TerHormonal <- read_dta("datasets/hersdata.dta")
# View(TerHormonal)
Variable Nombre Descripción Códigos/Valores
1 tchol1 Colesterol total al año de traramiento mg/dL
2 age Edad Años
3 raceth Raza 1=blanca, 2=afroamericana, 3=otra raza
4 BMI Índice de masa corporal kg/m^2
5 globrat Puntaje subjetivo de autoreporte de salud 1=pésimo estado de salur, 2=mal estado de salud, 3=regular estado de salud, 4=buen estado de salud, 5=estado de salud “ideal”
6 smoking Tabaquismo 0=No, 1=Si
7 physact Actividad física 1=sedentario, 2=actividad ocasional, 3=actividad regular, 4=activida frecuente, 5=actividad vigorosa
8 medcond Autorreporte enfermedad 0=No, 1=Si
9 diabetes Diabetes 0=No, 1=Si
10 tchol Colesterol total basal mg/dL
11 glucose Glucosa en sangre basal mg/dL
12 statins Uso de estatinas 0=No, 1=Si
13 LDL Niveles basales de LDL mg/dL
14 LDL1 Niveles de LDL al año mg/dL
15 HDL Niveles basales de HDL mg/dL
16 HDL1 Niveles de HDL al año mg/dL
17 TG Trigliceridos basales mg/dL
18 TG1 Trigliceridos al año de tratamiento mg/dL
19 SBP Presión arterial sistólica mmHg
20 DBP Presión arterial diastólica mmHg
# Esta práctica en el ejercicio real no se debe hacer.  Ver métodos de
# imputación de datos.

mydata <- na.omit(TerHormonal)

La base original tenía 2763 observaciones de 37 variables. Por fines académicos se excluyen registros con datos perdidos. Luego de eliminar observaciones con valores perdidos, quedando 2571 obsevaciones de 37 variables, es decir, se eliminaron 192 observacions (6.95% de los datos originales).

## Etiquetas de las variables

# Definiendo la raza como variable categorica
mydata$raceth <- factor(mydata$raceth, levels = c(1, 2, 3), labels = c("Blanca",
    "Afroamericana", "Otra raza"))

# Definiendo el puntaje subjetivo como categorica ordinal
mydata$globrat <- factor(mydata$globrat, levels = c(5, 4, 3, 2, 1), labels = c("estado de salud ideal",
    "buen estado de salud", "regular estado de salud", "mal estado de salud", "pésimo estado de salud"))

# Definiendo la variable tabaquismo como dicotomica
mydata$smoking <- factor(mydata$smoking, levels = c(0, 1), labels = c("No", "Si"))

# Definiendo actividad fisica como categorica ordinal
mydata$physact <- factor(mydata$physact, levels = c(5, 4, 3, 2, 1), labels = c("actividad vigorosa",
    "activida frecuente", "actividad regular", "actividad ocasional", "sedentario"))

# Definiendo la variable autoreporte de enfermedad como dicotomica
mydata$medcond <- factor(mydata$medcond, levels = c(0, 1), labels = c("No", "Si"))

# Definiendo la variable diabetes como dicotomica
mydata$diabetes <- factor(mydata$diabetes, levels = c(0, 1), labels = c("No", "Si"))

# Definiendo la variable uso de estatinas como dicotomica
mydata$statins <- factor(mydata$statins, levels = c(0, 1), labels = c("No", "Si"))

Análisis exploratorio de datos

A continuación se presenta un análisis exploratorio de las variables de interes:

library(summarytools)

# dfSummary(mydata,valid.col = FALSE, style = 'rmarkdown') # Permite hacer un
# análisis Para presentación en html

print(dfSummary(mydata, valid.col = FALSE), method = "render")

Data Frame Summary

mydata

Dimensions: 2571 x 37
Duplicates: 0
No Variable Label Stats / Values Freqs (% of Valid) Graph Missing
1 HT [haven_labelled, vctrs_vctr, double] random assignment to hormone therapy
Min : 0
Mean : 0.5
Max : 1
0:1303(50.7%)
1:1268(49.3%)
0 (0.0%)
2 age [numeric] age in years
Mean (sd) : 66.6 (6.6)
min ≤ med ≤ max:
44 ≤ 67 ≤ 79
IQR (CV) : 10 (0.1)
36 distinct values 0 (0.0%)
3 raceth [factor]
1. Blanca
2. Afroamericana
3. Otra raza
2299(89.4%)
184(7.2%)
88(3.4%)
0 (0.0%)
4 nonwhite [haven_labelled, vctrs_vctr, double] nonwhite race/ethnicity
Min : 0
Mean : 0.1
Max : 1
0:2299(89.4%)
1:272(10.6%)
0 (0.0%)
5 smoking [factor]
1. No
2. Si
2243(87.2%)
328(12.8%)
0 (0.0%)
6 drinkany [haven_labelled, vctrs_vctr, double] any current alcohol consumption
Min : 0
Mean : 0.4
Max : 1
0:1551(60.3%)
1:1020(39.7%)
0 (0.0%)
7 exercise [haven_labelled, vctrs_vctr, double] exercise at least 3 times per week
Min : 0
Mean : 0.4
Max : 1
0:1559(60.6%)
1:1012(39.4%)
0 (0.0%)
8 physact [factor]
1. actividad vigorosa
2. activida frecuente
3. actividad regular
4. actividad ocasional
5. sedentario
294(11.4%)
794(30.9%)
854(33.2%)
459(17.9%)
170(6.6%)
0 (0.0%)
9 globrat [factor]
1. estado de salud ideal
2. buen estado de salud
3. regular estado de salud
4. mal estado de salud
5. pésimo estado de salud
112(4.4%)
648(25.2%)
1229(47.8%)
536(20.8%)
46(1.8%)
0 (0.0%)
10 poorfair [haven_labelled, vctrs_vctr, double] poor/fair self-reported health
Min : 0
Mean : 0.2
Max : 1
0:1989(77.4%)
1:582(22.6%)
0 (0.0%)
11 medcond [factor]
1. No
2. Si
1624(63.2%)
947(36.8%)
0 (0.0%)
12 htnmeds [haven_labelled, vctrs_vctr, double] anti-hypertensive use
Min : 0
Mean : 0.8
Max : 1
0:464(18.0%)
1:2107(82.0%)
0 (0.0%)
13 statins [factor]
1. No
2. Si
1620(63.0%)
951(37.0%)
0 (0.0%)
14 diabetes [factor]
1. No
2. Si
1909(74.3%)
662(25.7%)
0 (0.0%)
15 dmpills [haven_labelled, vctrs_vctr, double] oral DM medication by self-report
Min : 0
Mean : 0.1
Max : 1
0:2325(90.4%)
1:246(9.6%)
0 (0.0%)
16 insulin [haven_labelled, vctrs_vctr, double] insulin use by self-report
Min : 0
Mean : 0.1
Max : 1
0:2327(90.5%)
1:244(9.5%)
0 (0.0%)
17 weight [numeric] weight (kg)
Mean (sd) : 72.7 (14.6)
min ≤ med ≤ max:
37.6 ≤ 70.9 ≤ 132
IQR (CV) : 19 (0.2)
541 distinct values 0 (0.0%)
18 BMI [numeric] BMI (kg/m^2)
Mean (sd) : 28.5 (5.4)
min ≤ med ≤ max:
15.5 ≤ 27.8 ≤ 54.1
IQR (CV) : 7 (0.2)
1400 distinct values 0 (0.0%)
19 waist [numeric] waist (cm)
Mean (sd) : 91.6 (13.5)
min ≤ med ≤ max:
56.9 ≤ 90.5 ≤ 170
IQR (CV) : 18 (0.1)
484 distinct values 0 (0.0%)
20 WHR [numeric] waist/hip ratio
Mean (sd) : 0.9 (0.1)
min ≤ med ≤ max:
0.6 ≤ 0.9 ≤ 1.2
IQR (CV) : 0.1 (0.1)
378 distinct values 0 (0.0%)
21 glucose [numeric] fasting glucose (mg/dl)
Mean (sd) : 111.7 (36.2)
min ≤ med ≤ max:
29 ≤ 99 ≤ 298
IQR (CV) : 23 (0.3)
198 distinct values 0 (0.0%)
22 weight1 [numeric] year 1 weight (kg)
Mean (sd) : 71.9 (14.9)
min ≤ med ≤ max:
37.7 ≤ 70.3 ≤ 142
IQR (CV) : 19.7 (0.2)
539 distinct values 0 (0.0%)
23 BMI1 [numeric] year 1 BMI (kg/m^2)
Mean (sd) : 28.3 (5.5)
min ≤ med ≤ max:
14.7 ≤ 27.5 ≤ 54
IQR (CV) : 7.2 (0.2)
1434 distinct values 0 (0.0%)
24 waist1 [numeric] year 1 waist (cm)
Mean (sd) : 91 (13.5)
min ≤ med ≤ max:
59 ≤ 90 ≤ 142
IQR (CV) : 19 (0.1)
488 distinct values 0 (0.0%)
25 WHR1 [numeric] year 1 waist/hip ratio
Mean (sd) : 0.9 (0.1)
min ≤ med ≤ max:
0.6 ≤ 0.9 ≤ 1.1
IQR (CV) : 0.1 (0.1)
366 distinct values 0 (0.0%)
26 glucose1 [numeric] year 1 fasting glucose (mg/dl)
Mean (sd) : 114.3 (44)
min ≤ med ≤ max:
42 ≤ 100 ≤ 440
IQR (CV) : 25 (0.4)
228 distinct values 0 (0.0%)
27 tchol [numeric] total cholesterol (mg/dl)
Mean (sd) : 228.2 (40.9)
min ≤ med ≤ max:
110 ≤ 224 ≤ 465
IQR (CV) : 51.5 (0.2)
222 distinct values 0 (0.0%)
28 LDL [numeric] LDL cholesterol (mg/dl)
Mean (sd) : 144.8 (37.8)
min ≤ med ≤ max:
36.8 ≤ 141 ≤ 393.4
IQR (CV) : 46.1 (0.3)
743 distinct values 0 (0.0%)
29 HDL [numeric] HDL cholesterol (mg/dl)
Mean (sd) : 50.3 (13.1)
min ≤ med ≤ max:
14 ≤ 49 ≤ 130
IQR (CV) : 16 (0.3)
84 distinct values 0 (0.0%)
30 TG [numeric] triglycerides (mg/dl)
Mean (sd) : 165.3 (63.1)
min ≤ med ≤ max:
31 ≤ 157 ≤ 476
IQR (CV) : 92 (0.4)
282 distinct values 0 (0.0%)
31 tchol1 [numeric] year 1 total cholesterol (mg/dl)
Mean (sd) : 219.2 (41.1)
min ≤ med ≤ max:
92 ≤ 214 ≤ 535
IQR (CV) : 49 (0.2)
229 distinct values 0 (0.0%)
32 LDL1 [numeric] year 1 LDL cholesterol (mg/dl)
Mean (sd) : 132.4 (39.1)
min ≤ med ≤ max:
-20 ≤ 128.8 ≤ 450.2
IQR (CV) : 47.6 (0.3)
751 distinct values 0 (0.0%)
33 HDL1 [numeric] year 1 HDL cholesterol (mg/dl)
Mean (sd) : 51.8 (13.9)
min ≤ med ≤ max:
14 ≤ 50 ≤ 124
IQR (CV) : 17 (0.3)
88 distinct values 0 (0.0%)
34 TG1 [numeric] year 1 triglycerides (mg/dl)
Mean (sd) : 175.2 (82.4)
min ≤ med ≤ max:
31 ≤ 157 ≤ 1010
IQR (CV) : 94 (0.5)
362 distinct values 0 (0.0%)
35 SBP [numeric] systolic blood pressure
Mean (sd) : 134.8 (19)
min ≤ med ≤ max:
83 ≤ 134 ≤ 224
IQR (CV) : 25 (0.1)
109 distinct values 0 (0.0%)
36 DBP [numeric] diastolic blood pressure
Mean (sd) : 73.1 (9.6)
min ≤ med ≤ max:
45 ≤ 72 ≤ 102
IQR (CV) : 13 (0.1)
58 distinct values 0 (0.0%)
37 age10 [numeric] age (per 10 years)
Mean (sd) : 6.7 (0.7)
min ≤ med ≤ max:
4.4 ≤ 6.7 ≤ 7.9
IQR (CV) : 1 (0.1)
36 distinct values 0 (0.0%)

Generated by summarytools 1.0.0 (R version 4.1.1)
2022-02-27

Correlaciones y gráficos de dispersión

library(PerformanceAnalytics)

# Seleccionar variables cuantitativas
mynumvar <- c("tchol1", "age", "BMI", "weight", "tchol", "glucose", "LDL", "HDL",
    "TG", "SBP", "DBP")

mydatanum <- mydata[, mynumvar]

chart.Correlation(mydatanum, histogram = FALSE, method = "pearson")

Según nuestro grafico de correlación hay alta correlación entre valor de colesterol total basal y colesterol total al año de tratamiento de 0,61. Desde el punto de vista medico no considero que haya una relación directa, incluso si se considera que es debido a la terapia hormonal, pues se espera que con la intervención incluso haya disminucion del valor basal.

El resto de los indices de correlación entre las demás variables y los niveles de colesterol a un año son menores de 0,6, lo que indica una correlacion moderada a lo sumo.

Matriz de correlaciones. Correlograma

library(corrplot)

corrmatriz <- as.matrix(round(cor(mydatanum), 2))

corrplot(corrmatriz, method = "number")

Gráficos de dispersión para el desenlace

Colesterol total al año e indice de masa corporal(BMI)

ggplot(mydata, aes(x = BMI, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS indice de masa corporal(BMI)") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que los datos son dispersos y que la pendiente es casi cercana a 0, indicando que la regresión lineal entre el BMI y los niveles de colesterol al año no es significativa, es decir el aumento en el BMI no afecta los niveles de colesterol al año.

Colesterol total al año y peso

ggplot(mydata, aes(x = weight, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS peso") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que los datos son dispersos y que la pendiente es casi cercana a 0, es decir el aumento en el peso no afecta los niveles de colesterol al año.

Colesterol total al año y edad en años

ggplot(mydata, aes(x = age, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS edad en años") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que los datos son dispersos y que la pendiente es casi cercana a 0, indicando que la regresión lineal entre la edad y los niveles de colesterol al año no es significativa, es decir el aumento en la edad no afecta los niveles de colesterol al año.

Colesterol total al año y colesterol basal

ggplot(mydata, aes(x = tchol, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS colesterol basal") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que hay una relación positiva en los niveles de colesterol total basal y los niveles de colesterol total al año con una pendiente positiva que indica que a mayor colesterol basal, mayor serán los niveles de colesterol al año, aunque en la dispersión de los datos observamos que esta relación descrita no es fuerte dado la distancia entre los datos y la recta. Además obervamos que existe una mayor concentración de los datos entre menos niveles de colesterol en ambas variable, con mayor dispersión de los mismos al final y con datos extremos que pueden ayudar a observar dicha dirección ya que como se dijo antes en el gráfico de correlación desde el punto de plausibilidad no es clara la relación descrita.

Colesterol total al año y actividad fisica

ggplot(mydata, aes(x = physact, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS actividad fisica") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

Los datos en esta gráfico de dispersión no muestra relación entre los niveles de actividad física y los niveles de colesterol total al año.

ggplot(mydata, aes(x = globrat, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS Puntaje de autorreporte de salud") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

Los datos en esta gráfico de dispersión no muestra relación entre el puntaje de autorreporte de salud y los niveles de colesterol total al año.

ggplot(mydata, aes(x = smoking, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS el tabaquismo") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

Los datos en esta gráfico de dispersión no muestra relación entre el tabquismo y los niveles de colesterol total al año.

ggplot(mydata, aes(x = medcond, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS el autorreporte de enfermedad") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

Los datos en esta gráfico de dispersión no muestra relación entre el autorreporte de enfermedad y el de colesterol total al año.

ggplot(mydata, aes(x = diabetes, y = tchol1)) + geom_point() + theme_classic() +
    ggtitle("Colesterol total al año VS la presencia de diabetes") + geom_smooth(method = "lm",
    se = FALSE, col = "red")

Los datos en esta gráfico de dispersión no muestra relación entre la presencia de diabetes y el de colesterol total al año.

ggplot(mydata, aes(x = glucose, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS niveles de glucosa") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que los datos son dispersos, con algunos datos extremos y que la pendiente es casi cercana a 0, indicando que la regresión lineal entre los niveles de glucosa y los niveles de colesterol al año no es significativa.

ggplot(mydata, aes(x = statins, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS uso de estatinas") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

Los datos en esta gráfico de dispersión no muestra relación entre el uso de estatinas y el de colesterol total al año. Aunque desde el punto de vista médico se debería observar una mayor concentración de colesterol total al año entre los aquellos que no usan estatinas.

ggplot(mydata, aes(x = LDL, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS niveles basales de LDL") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que hay una relación positiva en los niveles de colesterol LDL basal y los niveles de colesterol total al año con una pendiente positiva que indica que a mayor colesterol LDL basal, mayor serán los niveles de colesterol al año, similar a como ocurría en la gráfica que mencionamos de colesterol total basal y colesterol total al año, esto dado porque para el cálculo del colesterol total se requiere tener en cuenta el colesterol LDL. Por tanto se hace la misma anotación al respecto: desde el punto de plausibilidad no es clara la relación descrita.

ggplot(mydata, aes(x = LDL1, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS niveles de LDL al año") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

El en gráfico de dispersión observamos una relación lineal positiva y fuerte entre el grupo de datos, explicado porque para el cálculo de la variable de interes dependiente se requiere conocer el valor del LDL.

ggplot(mydata, aes(x = HDL, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS niveles de HDL basales") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que la relación lineal entre los niveles de HDL basales y los niveles de colesterol al año no es significativa.

ggplot(mydata, aes(x = HDL, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS niveles de HDL al año") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que la relación lineal entre los niveles de HDL al año y los niveles de colesterol al año no es significativa con una leve tendencia a una pendiente positiva.

ggplot(mydata, aes(x = TG, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS niveles de TG basales") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que la relación lineal entre los niveles de TG basales y los niveles de colesterol al año no es significativa.

ggplot(mydata, aes(x = TG1, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS niveles de TG al año") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que la relación lineal entre los niveles de TG al año y los niveles de colesterol al año no es significativa con una leve tendencia a una pendiente positiva pero explicada por datos extremos y una clara concentración mayor de datos.

ggplot(mydata, aes(x = SBP, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS la presión areterial sistólica") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que los datos son dispersos, con pendiente es casi cercana a 0, indicando que la relación lineal entre SBP y los niveles de colesterol al año no es significativa.

ggplot(mydata, aes(x = DBP, y = tchol1)) + geom_point() + theme_classic() + ggtitle("Colesterol total al año VS la presión areterial diastólica") +
    geom_smooth(method = "lm", se = FALSE, col = "red")

En el diagrama de dispersión podemos observar que los datos son dispersos, con pendiente es casi cercana a 0, indicando que la relación lineal entre DBP y los niveles de colesterol al año no es significativa.

Preguntas

1. Según los resultados anteriores, interprete y analice la relación entre colesterol total al año y los predictores cuantitativos.

A continuación se presentea una table de resumen en el que se indican la correlación y dispersión de las variables para su análisis. Considero que a la luz de los gráficos, análisis e índices descritos, en general hay una mala correlación lineal entre colesterol total al año y todos los predictores cuantitativos.

Tabla de resumen de correlación y dispersión variablescuantitativas

2. Las variables puntaje subjetivo del estado de salud (globrat) y actividad física (physact) son variables de naturaleza ordinal. ¿Qué coeficiente de correlación sugiere más adecuado para la relación con la variable dependiente?

Los coeficientes de correlación que se pueden implementar en este caso con variables policotómicas son dos: Coeficiente de correlación ETA. Coeficiente de correlación biserial por rangos.

Ambas se dan por regresión logística, el primero funciona con muchos niveles en variables ordinales, y el segundo busca convertir la variable policotómica en dicotómica reagrupando los grupos en dos variables nuevas. Por ejemplo, para la variable globrat, se crearía una categoría que agrupa los niveles 1, 2, 3 (mal o regular estado de salud) y otra que agrupa los niveles 4 y 5 (buen o muy buen estado de salud) y luego se calcula el coefieciente biserial por rangos.

Tabla ANOVA para el modelo lineal múltiple

modelo1 <- lm(tchol1 ~ age + BMI + weight + tchol + LDL + HDL + TG + SBP + DBP +
    glucose + globrat + physact + medcond + statins + raceth + smoking + diabetes,
    data = mydata)
summary(modelo1)

Call:
lm(formula = tchol1 ~ age + BMI + weight + tchol + LDL + HDL + 
    TG + SBP + DBP + glucose + globrat + physact + medcond + 
    statins + raceth + smoking + diabetes, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max 
-175.43  -20.12   -2.38   17.72  232.20 

Coefficients:
                                Estimate Std. Error t value Pr(>|t|)    
(Intercept)                     7.84e+01   1.13e+01    6.96  4.4e-12 ***
age                            -1.05e-01   1.09e-01   -0.96    0.335    
BMI                             5.83e-01   3.11e-01    1.87    0.061 .  
weight                         -2.58e-01   1.16e-01   -2.22    0.026 *  
tchol                           6.57e+04   1.68e+05    0.39    0.696    
LDL                            -6.57e+04   1.68e+05   -0.39    0.696    
HDL                            -6.57e+04   1.68e+05   -0.39    0.696    
TG                             -1.31e+04   3.36e+04   -0.39    0.696    
SBP                             9.42e-03   4.27e-02    0.22    0.825    
DBP                             1.35e-01   8.32e-02    1.62    0.105    
glucose                        -4.74e-02   2.53e-02   -1.87    0.062 .  
globratbuen estado de salud    -7.14e-01   3.36e+00   -0.21    0.832    
globratregular estado de salud -2.01e+00   3.29e+00   -0.61    0.542    
globratmal estado de salud     -4.14e+00   3.58e+00   -1.16    0.247    
globratpésimo estado de salud  -4.96e+00   6.04e+00   -0.82    0.412    
physactactivida frecuente       2.83e+00   2.24e+00    1.26    0.207    
physactactividad regular        3.24e+00   2.27e+00    1.43    0.154    
physactactividad ocasional      1.90e+00   2.59e+00    0.74    0.462    
physactsedentario               3.36e+00   3.39e+00    0.99    0.323    
medcondSi                       1.10e+00   1.40e+00    0.79    0.432    
statinsSi                       1.73e+00   1.38e+00    1.26    0.209    
racethAfroamericana             6.89e-01   2.60e+00    0.27    0.791    
racethOtra raza                 1.31e+00   3.64e+00    0.36    0.719    
smokingSi                       3.71e-01   2.03e+00    0.18    0.855    
diabetesSi                      4.18e+00   2.14e+00    1.96    0.050 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 32.6 on 2546 degrees of freedom
Multiple R-squared:  0.38,  Adjusted R-squared:  0.374 
F-statistic:   65 on 24 and 2546 DF,  p-value: <2e-16

3. El coeficiente de determinación. ¿Cuando seleccionar Multiple R-square o Adjusted R-squared? y ¿Por qué?. Interprete el coeficiente de determinación del modelo anterior.

El coeficiente de correlación (R2) es el porcentaje de variación de la variable dependiente (niveles colesterol total al año) que se explica por el modelo. Por lo general, mientras mayor sea el R2, mejor será el ajuste del modelo a sus datos y los valores van entre 0 y 1. El R2 ajustado es el porcentaje de variación en la variable de respuesta que es explicado por su relación con una o más variables predictoras, ajustado para el número de predictores en el modelo, es decir, en este el el porcentaje de variación no cambiará independientemente de que se incluyan mas variables al modelo, dando una clara idea del modelo.

Para nuestro ejemplo, el modelo explica el 38% de los resultados predichos en la variable dependiente (niveles colesterol total al año). Lo que lo convierte en un modelo con efectividad mala o regular.

4a. Interprete 1-R^2

1 - 38% = 62%.

Indica el porcentaje de la variación de la variable de respuesta no explicado por el modelo. Esto incluye tanto variables no conocidas o no incluidas como el azar. En este caso especifico pudieran incluirse variables relacionadas con la dieta que puedan en un momento dado explicar los niveles de colesterol total al año.

1-R^2 se interpreta como el porcentaje de variación en la variable de respuesta que NO es explicado por su relación con una o más variables predictoras, ajustado para el número de predictores en el modelo.

4b. ¿Qué otras variables incluiría en el modelo desde su conocimiento que explique la variabilidad del colesterol total al año?

# Tabla de los coeficientes del modelo
modelo1 <- lm(tchol1 ~ age + BMI + weight + tchol + LDL + HDL + TG + SBP + DBP +
    glucose + globrat + physact + medcond + statins + raceth + smoking + diabetes,
    data = mydata)
# summary(modelo1)

library(sjPlot)  #Para la presentación de la tabla con los coeficientes
tab_model(modelo1, show.se = TRUE, show.stat = TRUE, dv.labels = c("Modelo 1"))
  Modelo 1
Predictors Estimates std. Error CI Statistic p
(Intercept) 78.35 11.26 56.27 – 100.44 6.96 <0.001
age in years -0.11 0.11 -0.32 – 0.11 -0.96 0.335
BMI (kg/m^2) 0.58 0.31 -0.03 – 1.19 1.87 0.061
weight (kg) -0.26 0.12 -0.49 – -0.03 -2.22 0.026
total cholesterol (mg/dl) 65709.58 167884.34 -263494.18 – 394913.34 0.39 0.696
LDL cholesterol (mg/dl) -65708.97 167884.34 -394912.73 – 263494.79 -0.39 0.696
HDL cholesterol (mg/dl) -65708.92 167884.34 -394912.68 – 263494.85 -0.39 0.696
triglycerides (mg/dl) -13141.80 33576.87 -78982.55 – 52698.95 -0.39 0.696
systolic blood pressure 0.01 0.04 -0.07 – 0.09 0.22 0.825
diastolic blood pressure 0.14 0.08 -0.03 – 0.30 1.62 0.105
fasting glucose (mg/dl) -0.05 0.03 -0.10 – 0.00 -1.87 0.062
globrat: buen estado de
salud
-0.71 3.36 -7.30 – 5.87 -0.21 0.832
globrat: regular estado
de salud
-2.01 3.29 -8.46 – 4.45 -0.61 0.542
globrat: mal estado de
salud
-4.14 3.58 -11.15 – 2.87 -1.16 0.247
globrat: pésimo estado de
salud
-4.96 6.04 -16.81 – 6.90 -0.82 0.412
physact: activida
frecuente
2.83 2.24 -1.57 – 7.23 1.26 0.207
physact: actividad
regular
3.24 2.27 -1.22 – 7.69 1.43 0.154
physact: actividad
ocasional
1.90 2.59 -3.17 – 6.98 0.74 0.462
physact: sedentario 3.36 3.39 -3.30 – 10.01 0.99 0.323
medcond: Si 1.10 1.40 -1.64 – 3.84 0.79 0.432
statins: Si 1.73 1.38 -0.97 – 4.44 1.26 0.209
raceth: Afroamericana 0.69 2.60 -4.41 – 5.79 0.27 0.791
raceth: Otra raza 1.31 3.64 -5.83 – 8.44 0.36 0.719
smoking: Si 0.37 2.03 -3.60 – 4.34 0.18 0.855
diabetes: Si 4.18 2.14 -0.01 – 8.37 1.96 0.050
Observations 2571
R2 / R2 adjusted 0.380 / 0.374

Evaluando los coeficientes para el modelo lineal muúltiple: Intercepto es significativo BMI (valor p=0.061 muy cercano a 0.05) Weight (valor p=0.026) Glucose (valor p=0.062 cercano a 0.05) Diabetes (valor p=0.050)

Evaluando las Correlaciones con la respuesta: BMI Correlación muy baja weight Correlación muy baja glucose Correlación muy baja

Evaluando la multicolinealidad: corr entre regresores BMI vs weight Corr muy alta, en el modelo final incluyo únicamente una de ellas preferiblemente la que tenga mayor corr con la respuesta. BMI con glucose Corr baja, se pueden incluir ambas regresoras. Weight con glucose Corr baja, se pueden incluir ambas regresoras.

Modelo seleccionado Basados únicamente en los tres ítems anteriores (Coeficientes, Correlación y Colinealidad) el modelo sería:

Col_1año = 78.4 + 0.583xBMI - 0.0474xglucose + 4.18xdiabetes_Si

Plausibilidad médica para incluir otras variables

tchol porque el tchol1 depende puede influencierse directamente por los niveles baselaes de tchol y verse reflejado en el valor del mismo a un año. LDL porque se requiere el mismo para el cálculo del colesterol total. HDL y tchol1 tienen una relación inversa en cuanto al impacto en la salud cardiovascular. physactividad porque impacta positivamente los niveles de colesterol en la sangre. statins porque este medicamento regula tchol1..

tchol1 vs tchol Corr alta 0.61 tchol1 vs LDL Corr alta 0.57 tchol1 vs HDL Corr muy baja

tchol vs LDL Corr muy alta tchol vs HDL Corr muy baja LDL vs HDL Corr muy baja

En conclusión, por plausabilidad médica incluyo tchol nada más, y el modelo queda:

Col_1año = 78.4 + 0.583xBMI - 0.0474xglucose + 4.18xdiabetes_Si + 65700xtchol

Sin embargo, notese que el coeficiente de tchol “dispara” el valor de tchol1 a valores muy lejanos de su rango y su IC, por lo que es mejor no incluir esa variable.

Interpretación de los coeficientes del modelo

tab_model(modelo1, show.se = TRUE, show.stat = TRUE, dv.labels = c("Modelo 1"))
  Modelo 1
Predictors Estimates std. Error CI Statistic p
(Intercept) 78.35 11.26 56.27 – 100.44 6.96 <0.001
age in years -0.11 0.11 -0.32 – 0.11 -0.96 0.335
BMI (kg/m^2) 0.58 0.31 -0.03 – 1.19 1.87 0.061
weight (kg) -0.26 0.12 -0.49 – -0.03 -2.22 0.026
total cholesterol (mg/dl) 65709.58 167884.34 -263494.18 – 394913.34 0.39 0.696
LDL cholesterol (mg/dl) -65708.97 167884.34 -394912.73 – 263494.79 -0.39 0.696
HDL cholesterol (mg/dl) -65708.92 167884.34 -394912.68 – 263494.85 -0.39 0.696
triglycerides (mg/dl) -13141.80 33576.87 -78982.55 – 52698.95 -0.39 0.696
systolic blood pressure 0.01 0.04 -0.07 – 0.09 0.22 0.825
diastolic blood pressure 0.14 0.08 -0.03 – 0.30 1.62 0.105
fasting glucose (mg/dl) -0.05 0.03 -0.10 – 0.00 -1.87 0.062
globrat: buen estado de
salud
-0.71 3.36 -7.30 – 5.87 -0.21 0.832
globrat: regular estado
de salud
-2.01 3.29 -8.46 – 4.45 -0.61 0.542
globrat: mal estado de
salud
-4.14 3.58 -11.15 – 2.87 -1.16 0.247
globrat: pésimo estado de
salud
-4.96 6.04 -16.81 – 6.90 -0.82 0.412
physact: activida
frecuente
2.83 2.24 -1.57 – 7.23 1.26 0.207
physact: actividad
regular
3.24 2.27 -1.22 – 7.69 1.43 0.154
physact: actividad
ocasional
1.90 2.59 -3.17 – 6.98 0.74 0.462
physact: sedentario 3.36 3.39 -3.30 – 10.01 0.99 0.323
medcond: Si 1.10 1.40 -1.64 – 3.84 0.79 0.432
statins: Si 1.73 1.38 -0.97 – 4.44 1.26 0.209
raceth: Afroamericana 0.69 2.60 -4.41 – 5.79 0.27 0.791
raceth: Otra raza 1.31 3.64 -5.83 – 8.44 0.36 0.719
smoking: Si 0.37 2.03 -3.60 – 4.34 0.18 0.855
diabetes: Si 4.18 2.14 -0.01 – 8.37 1.96 0.050
Observations 2571
R2 / R2 adjusted 0.380 / 0.374

5. Interprete el coeficiente de regresión de: edad, IMC (BMI), peso (weight), auto reporte de salud (self-reported) y colesterol basal.

Interpretación:

** Edad: Por un incremento en un año de edad de los individuos, los niveles de colesterol total al año disminuye en 0.11mg/dL, despues de ajustar por IMC, peso, niveles de colesterol total basal, niveles de LDL basal, niveles de HDL basal, niveles de TG total basal, niveles de glucosa, presión arterial sistólica, presión arterial diastólica, puntaje subjetivo de autorreporte de salud, actividad física, consumo de estatinas, raza, fumar y diabetes.

** IMC (BMI): Por un incremento en 1kg/m^2 de los individuos, los niveles de colesterol total al año aumentan en 0.58mg/dL, despues de ajustar por edad, peso, niveles de colesterol total basal, niveles de LDL basal, niveles de HDL basal, niveles de TG total basal, niveles de glucosa, presión arterial sistólica, presión arterial diastólica, puntaje subjetivo de autorreporte de salud, actividad física, consumo de estatinas, raza, fumar y diabetes.

** Auto reporte de salud (self-reported): Siempre que los individuos reportan buen estado de salud, los niveles de colesterol total al año disminuye en 0.71mg/dL, si reportan regular estado de salud disminuye en 2.01mg/dL, si reportan mal estado de salud disminuye en 4.14mg/dL, si reportan estado de salud pesimo disminuye en 4.96mg/dL. Todo esto luego de ajustar por edad, IMC, peso, niveles de colesterol total basal, niveles de LDL basal, niveles de HDL basal, niveles de TG total basal, niveles de glucosa, presión arterial sistólica, presión arterial diastólica, actividad física, consumo de estatinas, raza, fumar y diabetes.

** Colesterol basal: Por un incremento en 1mg/dL en los niveles de colesterol total basal de los individuos, los niveles de colesterol total al año aumentan en 65709.58mg/dL, IMC, peso, niveles de LDL basal, niveles de HDL basal, niveles de TG total basal, niveles de glucosa, presión arterial sistólica, presión arterial diastólica, puntaje subjetivo de autorreporte de salud, actividad física, consumo de estatinas, raza, fumar y diabetes. Sin embargo estos valores están por fuera de los rangos de los datos de la base de datos de los pacientes, además, no son plaucibes biológicamente.

6. Desde su punto de vista analice la precisión de los intervalos de confianza.

En cuanto a los intervalos de confianza de los coeficientes de regresión para las variables descritas, todos cruzan el cero. Para la edad y para el IMC, los intervalos son pequeños y su cambio no tiene mucha importancia desde el punto de vista clínico en cuanto a la variable de interes (tchol1), mientras que para el estado de salud autoreportado los intervalos de confianza son amplios mayor cambio en los niveles de colesterol total al año. Para el colesterol, como se mencionó anteriormente se presentarón datos irreales y con intervalos muy amplios por fuera de los rangos minimos y maximos presentados en la base de datos original.

7. Desde el enfoque de la significación estadística (cada vez más en desuso: ver discusión del valor p), asumiendo un alfa=0.05, ¿Qué coeficientes son estadísticamente significativos? (Beta diferente de Cero como hipótesis alterna). En sus palabras: ¿Porqué este enfoque puede ser limitado?

Los coeficientes estadisticamente significativos son aquellos con p<0.05, en nuestro modelo el intercepto, la variable peso (weight) y la presencia de diabetes.

El enfoque del valor la significancia estadística y el valor de “p”, es un enfoque reduccionista porque no se deben basar las conclusiones en estos resultados solo por tener un valor de p<0.05, o no se puede pretender que no hay efecto porque no se presentó dicho resultado. Siempre habrá espacio para la incertidumbre y el azar en la investigación.

Análisis del problema de la confusión

8. Compare la magnitud y dirección del efecto (coeficiente beta asociado al ejercicio) crudo o sin ajustar con el efecto ajustado para la relación ejercicio y colesterol total al año. ¿Qué otras variables potenciales de confusión incluiría en el modelo para la relación ejercicio y colesterol al año?

Variable confusora: BMI

Comenzamos calculando el modelo crudo de las variables ejercicio y colesterol total al año. Luego, si suponemos que la variable de confusión es BMI, calculemos el modelo ajustado por indice de masa corporal del modelo anterior:

# Modelo ejercicio y colesterol total al año (crudo)
modeloCrudo <- lm(tchol1 ~ exercise, data = mydata)

# Modelo ejercicio y colesterol total al año (crudo) ajustado por BMI
modeloAjustadoBMI <- lm(tchol1 ~ exercise + BMI, data = mydata)

# Mostrando ambas tablas
tab_model(modeloCrudo, modeloAjustadoBMI, show.se = TRUE, show.stat = TRUE, dv.labels = c("Modelo Crudo tchol1 vs exercise",
    "Modelo ajustado por BMI"))
  Modelo Crudo tchol1 vs exercise Modelo ajustado por BMI
Predictors Estimates std. Error CI Statistic p Estimates std. Error CI Statistic p
(Intercept) 219.68 1.04 217.63 – 221.72 210.78 <0.001 214.45 4.55 205.52 – 223.37 47.12 <0.001
exercise at least 3 times
per week
-1.15 1.66 -4.41 – 2.11 -0.69 0.489 -0.81 1.69 -4.11 – 2.50 -0.48 0.631
BMI (kg/m^2) 0.18 0.15 -0.12 – 0.47 1.18 0.238
Observations 2571 2571
R2 / R2 adjusted 0.000 / -0.000 0.001 / -0.000

Tanto en el modelo crudo como en el ajustado la pendiente (coeficiente de exercise) no dió significativa, y además el coeficiente R2/R2 ajustado dió muy bajo, por lo que el modelo lineal entre niveles de coresterol al año y frecuencia de ejercicio no es significativo. BMI en el modelo ajustado provocó un aumento en la pendiente de -1.15 en el modelo crudo a -0.81 en el modelo ajustado, aunque ambos valores muy cercanos a cero. En general La inclusión de BMI en el modelo no afecta para nada la significancia de del modelo lineal.

Variable confusora: physact

Comenzamos calculando el modelo crudo de las variables ejercicio y colesterol total al año. Luego, si suponemos que la variable de confusión es physactividad, calculemos el modelo ajustado por actividad física corporal del modelo anterior:

# Modelo ejercicio y colesterol total al año (crudo)
modeloCrudo <- lm(tchol1 ~ exercise, data = mydata)

# Modelo ejercicio y colesterol total al año (crudo) ajustado por physact
modeloAjustadophisactividad <- lm(tchol1 ~ exercise + physact, data = mydata)

# Mostrando ambas tablas
tab_model(modeloCrudo, modeloAjustadophisactividad, show.se = TRUE, show.stat = TRUE,
    dv.labels = c("Modelo Crudo tchol1 vs exercise", "Modelo ajustado por physact"))
  Modelo Crudo tchol1 vs exercise Modelo ajustado por physact
Predictors Estimates std. Error CI Statistic p Estimates std. Error CI Statistic p
(Intercept) 219.68 1.04 217.63 – 221.72 210.78 <0.001 217.21 2.63 212.05 – 222.36 82.60 <0.001
exercise at least 3 times
per week
-1.15 1.66 -4.41 – 2.11 -0.69 0.489 -1.21 1.74 -4.63 – 2.21 -0.69 0.488
physact: activida
frecuente
3.40 2.82 -2.13 – 8.92 1.21 0.228
physact: actividad
regular
3.27 2.82 -2.25 – 8.80 1.16 0.245
physact: actividad
ocasional
1.37 3.16 -4.83 – 7.56 0.43 0.665
physact: sedentario 1.69 4.05 -6.26 – 9.64 0.42 0.677
Observations 2571 2571
R2 / R2 adjusted 0.000 / -0.000 0.001 / -0.001

Tanto en el modelo crudo como en el ajustado la pendiente (coeficiente de exercise) no dió significativa y el coeficiente R2/R2 ajustado dió muy bajo, por lo que el modelo lineal entre niveles de coresterol al año y frecuencia de ejercicio no es significativo. Con la inclusión de la variables actividad física no hay cambios en la pendiente y no afecta el modelo lineal.

Variable confusora: drinkany

# Modelo ejercicio y colesterol total al año (crudo)
modeloCrudo <- lm(tchol1 ~ exercise, data = mydata)

# Modelo ejercicio y colesterol total al año (crudo) ajustado por consumo de
# alcohol
modeloAjustadoAlcohol <- lm(tchol1 ~ exercise + drinkany, data = mydata)

# Mostrando ambas tablas
tab_model(modeloCrudo, modeloAjustadoAlcohol, show.se = TRUE, show.stat = TRUE, dv.labels = c("Modelo Crudo tchol1 vs exercise",
    "Modelo ajustado por consumo de alcohol"))
  Modelo Crudo tchol1 vs exercise Modelo ajustado por consumo de alcohol
Predictors Estimates std. Error CI Statistic p Estimates std. Error CI Statistic p
(Intercept) 219.68 1.04 217.63 – 221.72 210.78 <0.001 219.29 1.22 216.89 – 221.68 179.53 <0.001
exercise at least 3 times
per week
-1.15 1.66 -4.41 – 2.11 -0.69 0.489 -1.18 1.66 -4.44 – 2.08 -0.71 0.477
any current alcohol
consumption
1.01 1.66 -2.24 – 4.27 0.61 0.542
Observations 2571 2571
R2 / R2 adjusted 0.000 / -0.000 0.000 / -0.000

Tanto en el modelo crudo como en el ajustado la pendiente (coeficiente de exercise) no dió significativa y el coeficiente R2/R2 ajustado dió muy bajo, por lo que el modelo lineal entre niveles de coresterol al año y frecuencia de ejercicio no es significativo. Con la inclusión de la variables actividad física no hay camibos en la pendiente y no afecta el modelo lineal.

9. Plantee un ejemplo en su área de interés en el que se pueda usar un modelo de regresión lineal con fines explicativos. Indique cuál sería la variable dependiente, la variable independiente y las potenciales variables de confusión.

Se plantea un estudio para evaluar el impacto de la duración del tratamiento de hiperplasia prostática benigna en la función sexual (Medido en escalas función sexual). La variable dependiente será función sexual y dentro de las variables independientes está la residencia, el nivel de estudio, estrato socio económico, situación laboral, tratamiento, duración de tratamiento, comorbilidades, consumo de medicamentos, estado civil, edad. Dentro de las variables de confusión podemos encontrar comorbilidades que se presenten en pacientes mayores de 50 años o el consumo de medicamentos que pueden generar disfunción erectil.

10. Plantee un ejemplo en su área de interés en el que se pueda usar un modelo de regresión lineal con fines predictivos Indique cuál sería la variable dependiente, las variables predictoras.

Se plantea un estudio en infecciones urinarias en pacientes postoperatorios de cirugía urológica. Se busca predecir el riesgo de sepsis en el postoperatorio. La variable dependiente es el porcentaje de infecciones que se presentan luego de cirugía urológica y dentro de las variables predictoras se incluirá la presencia de infecciones urologicas previas a cirugía, la bacteriuria asintomática previa a cirugía, la presencia de sonda o cateter en el paciente, la edad, la obesidad, los días de hospitalización, el tipo de cirugía, el tiempo quirúrgico.