Proyecto 2 - Estadística avanzada.
Gabriel Missael Barco. NUA: 427071 | Laura del Carmen Cabal Paramo. NUA: 39158
21 de Octubre de 2021
Introducción.
Para la realización de este proyecto, se pretende reproducir de manera parcial (de acuerdo a los datos disponibles) los resultados del artículo:
Gallagher R, McKinley S, Dracup K. Predictors of women’s attendance at cardiac rehabilitation programs. Prog Cardiovasc Nurs. 2003 Summer;18(3):121-6. doi: 10.1111/j.0889-7204.2003.02129.x. PMID: 12893973.
De este artículo, se realizó un estudio descriptivo para identificar los factores que influyen en la asistencia de las mujeres a los programas de rehabilitación cardíaca y la adherencia de las mujeres a la modificación de los factores de riesgo después de un evento cardíaco.
Metodología seguida en artículo.
El estudio utilizó un diseño descriptivo con medidas al alta hospitalaria y de 12 semanas después del alta. Una muestra de conveniencia de pacientes mujeres (N = 196) fue reclutado de cuatro hospitales metropolitanos en Sydney, Australia. Los programas de rehabilitación cardiaca incluyeron sesiones de ejercicio supervisadas, prescripciones para hacer ejercicio en el hogar y componentes de educación sobre modificación de factores de riesgo. Los programas variaron en duración de 6 a 12 semanas Las mujeres rechazaron la inscripción debido a síntomas (n = 6), falta de interés (n = 5) y falta de voluntad para firmar un formulario de consentimiento (n = 4). La mayoría de las mujeres (n = 183) completaron las 12 semanas de estudio. Las mujeres que no completaron fueron 13, ya que habían muerto (n = 5), retirados (n = 4), o se perdieron durante el seguimiento (n = 4). Para la recolección de datos, se les llamó a estas mujeres y se les preguntaron si habían completado los programas y el por qué no si es que no lo habían hecho.
Resultados del artículo.
Las mujeres (N = 196) ingresadas en el hospital por un evento cardíaco se les dio seguimiento a las 12 semanas después del alta. Solo el 64% (n = 112) había sido referido a programas de rehabilitación cardíaca. A las 12 semanas posteriores al alta, solo el 32% de la muestra total (n = 57) asistía a programas y el 12% de la muestra total (n = 21) había abandonado antes de la finalización. Las probabilidades de que una mujer asistiera a rehabilitación cardíaca disminuyeron por el diagnóstico de infarto de miocardio, la falta de empleo, <55 años o> 70 años, y experimentar un evento personal estresante durante el seguimiento. Es probable que las mujeres se adhieran a las pautas de modificación del tabaquismo, la medicación y el estrés, pero es poco probable que se adhieran a las pautas de modificación de la dieta y el ejercicio.
Análisis de los datos en el artículo.
Se calculó un tamaño de muestra de 180 para el análisis de regresión logística, basado en un mínimo de 10 pacientes por celda para las 10 variables y variables ficticias (dummy variables) adicionales ingresadas. Se realizó un análisis de regresión logística sobre la asistencia utilizando las 10 variables: edad, educación, tipo de hospital, trabajo, arreglos de vivienda, diagnóstico, control percibido, reingreso, complicación, y evento estresante personal durante el seguimiento. Los resultados de estos análisis se presentan como odds ratios con nivel de confianza del 95%. Las diferencias en la modificación de los factores de riesgo a lo largo del tiempo se analizaron utilizando análisis de varianza de medidas repetidas. Dado este analisis presentado, vamos a replicar la regresión logística y odds ratio para este proyecto.
Reproducción parcial de resultados.
Regresión logistica sobre la edad.
# Cargamos los datos
df <- read.csv('dataset1.csv')
# Visualizamos estructura del dataset
head(df, n = 10)
# Libreria para graficar
library(ggplot2)
ggplot(df, aes(x= AGE, y=ATT, color=ATT)) +
geom_point(size=3)
Realizamos la regresión logistica:
model <- glm(ATT ~ AGE, family = binomial("logit"), data = df)
summary(model)
Call:
glm(formula = ATT ~ AGE, family = binomial("logit"), data = df)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.4345 -0.8975 -0.7646 1.3674 1.7589
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.87467 0.98088 1.911 0.05598 .
AGE -0.03788 0.01462 -2.590 0.00959 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 236.48 on 183 degrees of freedom
Residual deviance: 229.52 on 182 degrees of freedom
AIC: 233.52
Number of Fisher Scoring iterations: 4Por lo tanto, tenemos que la ecuación de regresión es:
\[\ln { \frac {p}{1-p} } = -0.03788x + 1.8746\]
Y la ecuación para estimar \(p\) es:
\[p(x) = \frac{e^{-0.03788x + 1.8746}}{1+e^{-0.03788x + 1.8746}}\] Podemos graficar esta función junto con los datos:
# Creamos arreglo para gráfica
x = c(0:100)
# Función de regresión
p_log <- function(x){
p <- (exp(-0.03788*x + 1.8746))/(1+exp(-0.03788*x + 1.8746))
}
# Probabilidades (odds)
y <- p_log(x)
# Graficamos la regresión
plot(x, y, type = 'l',
main = 'Logistic regresión and data',
col = 'blue', xlab = 'AGE',
ylab = 'p', ylim = c(0, 1))
# Graficamos los datos
lines(df$AGE, df$ATT, type = 'p', col = 'red')
# Etiquetas de líneas.
legend("center", c("Logistic regression", "Data"), fill=c("blue", 'red'))
# Agregamos una cuadrícula.
grid()
Podemos hacer una prueba de hipótesis de \(H_0: \beta_1 = 0\), esto es, que la pendiente de la regresión logistica es cero. Para ello, del summary anterior tenemos que el estadístico de prueba es \(z = -2.59\). Consideramos una significancia de \(\alpha = 0.01\), y el valor crítico:
qnorm(0.005)[1] -2.575829qnorm(1-0.005)[1] 2.575829Por lo tanto, dado que el estadístico de prueba es menor al valor crítico inferior, rechazamos la hipótesis nula \(H_0\) y concluimos que el modelo es una buena manera de explicar los datos.
El odds ratio en este caso:
exp(-0.03788)[1] 0.9628285Esto es, hay un \(3.8\%\) mas de probabilidad de que una persona no asista a las sesiones.
Regresión logistica sobre el indice de ansiedad y depresión.
# Cargamos los datos
df <- read.csv('dataset2.csv')
# Visualizamos estructura del dataset
head(df, n = 10)
# Libreria para graficar
library(ggplot2)
ggplot(df, aes(x= INDEX, y=PART, color=PART)) +
geom_point(size=3)
Realizamos la regresión logistica:
model <- glm(PART ~ INDEX, family = binomial("logit"), data = df)
summary(model)
Call:
glm(formula = PART ~ INDEX, family = binomial("logit"), data = df)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.2774 -0.9236 -0.7704 1.2777 1.9591
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.76027 0.45831 -3.841 0.000123 ***
INDEX 0.06641 0.02516 2.639 0.008308 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 236.48 on 183 degrees of freedom
Residual deviance: 229.14 on 182 degrees of freedom
AIC: 233.14
Number of Fisher Scoring iterations: 4Por lo tanto, tenemos que la ecuación de regresión es:
\[\ln { \frac {p}{1-p} } = 0.06641x -1.76027\]
Y la ecuación para estimar \(p\) es:
\[p(x) = \frac{e^{0.06641x -1.76027}}{1+e^{0.06641x -1.76027}}\] Podemos graficar esta función junto con los datos:
# Creamos arreglo para gráfica
x = c(-20:60)
# Función de regresión
p_log <- function(x){
p <- (exp(0.06641*x -1.76027))/(1+exp(0.06641*x -1.76027))
}
# Probabilidades (odds)
y <- p_log(x)
# Graficamos la regresión
plot(x, y, type = 'l',
main = 'Logistic regresión and data',
col = 'blue', xlab = 'INDEX',
ylab = 'p', ylim = c(0, 1))
# Graficamos los datos
lines(df$INDEX, df$PART, type = 'p', col = 'red')
# Etiquetas de líneas.
legend("center", c("Logistic regression", "Data"), fill=c("blue", 'red'))
# Agregamos una cuadrícula.
grid()
Podemos hacer una prueba de hipótesis de \(H_0: \beta_2= 0\), esto es, que la pendiente de la regresión logistica es cero. Para ello, del summary anterior tenemos que el estadístico de prueba es \(z = 2.639\). Consideramos una significancia de \(\alpha = 0.01\), y el valor crítico:
qnorm(0.005)[1] -2.575829qnorm(1-0.005)[1] 2.575829Por lo tanto, dado que el estadístico de prueba es menor al valor crítico inferior, rechazamos la hipótesis nula \(H_0\) y concluimos que el modelo es una buena manera de explicar los datos.
Podemos obtener también explícitamente la probabilidad de participar en el estudio dado el indice de estrés medido:
INDEX <- df$INDEX
P <- p_log(INDEX)
REAL <- df$PART
probabilities <- data.frame(INDEX, P, REAL)
head(probabilities, n = 5)
tail(probabilities, n = 5)El odds ratio en este caso:
exp(0.06641)[1] 1.068665Esto es, hay un \(6\%\) mas de probabilidad de que una persona asista a las sesiones.
Referencias.
- Gallagher R, McKinley S, Dracup K. Predictors of women’s attendance at cardiac rehabilitation programs. Prog Cardiovasc Nurs. 2003 Summer;18(3):121-6. doi: 10.1111/j.0889-7204.2003.02129.x. PMID: 12893973.
- Wayne W. Daniel y Chad L. Cross, Biostatistics: A Foundation for Analysis in the Health Sciences, 10° edition, Wiley.
Anexos.
Diapositivas de presentación de proyecto.
Se pueden encontrar en este link.