Primero vamos a abrir las librerías necesarias

library(readxl)  ## para abrir el archivo de excel
library(dplyr) ## para manipulación de las variables, ej revisar datos faltantes
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(psych) ## para realizar análisis descriptivo
library(Hmisc)# para otorgarle labels a las variables
## 
## Attaching package: 'Hmisc'
## The following object is masked from 'package:psych':
## 
##     describe
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(compareGroups) # para realizar tabla1
library(epiR) # para medidas de efecto
## Loading required package: survival
## Package epiR 2.0.78 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyses
## 
library(descr) # para tablas cruzadas
library(ggplot2) # para hacer gráficos
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
library(car) #para test de Levene
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:psych':
## 
##     logit
## The following object is masked from 'package:dplyr':
## 
##     recode

Ahora vamos a importar la base de datos con el nombre “BaseFram”, y vamos a abrirla para visualizarla

BaseFram <- read_excel("C:/Users/maria/OneDrive/Escritorio/iecs/TAKEHOME bio/Framingham (1).xlsx")
View(BaseFram)

Vamos a empezar a filtrar la base de datos para limpiarla. Primero vamos a seleccionar solo las variables que vamos a usar para el trabajo, y vamos a nombrar a esta base de datos filtrada “BaseTH”

BaseTH<- BaseFram[, c(2,4,5,7,8,9,10,14,17,18,21,22,23,25,26)]
print(BaseTH)
## # A tibble: 4,434 × 15
##      SEX   AGE SYSBP CURSMOKE CIGPDAY   BMI DIABETES PREVCHD PREVSTRK PREVHYP
##    <dbl> <dbl> <dbl>    <dbl>   <dbl> <dbl>    <dbl>   <dbl>    <dbl>   <dbl>
##  1     1    39  106         0       0  27.0        0       0        0       0
##  2     0    46  121         0       0  28.7        0       0        0       0
##  3     1    48  128.        1      20  25.3        0       0        0       0
##  4     0    61  150         1      30  28.6        0       0        0       1
##  5     0    46  130         1      23  23.1        0       0        0       0
##  6     0    43  180         0       0  30.3        0       0        0       1
##  7     0    63  138         0       0  33.1        0       0        0       0
##  8     0    45  100         1      20  21.7        0       0        0       0
##  9     1    52  142.        0       0  26.4        0       0        0       1
## 10     1    43  162         1      30  23.6        0       0        0       1
## # ℹ 4,424 more rows
## # ℹ 5 more variables: DEATH <dbl>, ANGINA <dbl>, HOSPMI <dbl>, ANYCHD <dbl>,
## #   STROKE <dbl>

A todas las variables con datos categóricos, vamos a convertirlas en factor

BaseTH$SEX<-as.factor(BaseTH$SEX)
BaseTH$CURSMOKE<-as.factor(BaseTH$CURSMOKE)
BaseTH$DIABETES<-as.factor(BaseTH$DIABETES)
BaseTH$PREVCHD<-as.factor(BaseTH$PREVCHD)
BaseTH$PREVHYP<-as.factor(BaseTH$PREVHYP)
BaseTH$PREVSTRK<-as.factor(BaseTH$PREVSTRK)
BaseTH$DEATH<-as.factor(BaseTH$DEATH)
BaseTH$ANGINA<-as.factor(BaseTH$ANGINA)
BaseTH$HOSPMI<-as.factor(BaseTH$HOSPMI)
BaseTH$ANYCHD<-as.factor(BaseTH$ANYCHD)
BaseTH$STROKE<-as.factor(BaseTH$STROKE)

Vamos a revisar si existen datos faltantes en las variables, y luego vamos a eliminarlos. Para esto, primero creamos un objeto llamado “observaciones_faltantes”, y luego las visualizo en la base de datos para identificar a qué variables pertenecen. Finalmente, filtro mi base eliminando las filas con datos faltantes.

observaciones_faltantes <- which(!complete.cases(BaseTH))
print (observaciones_faltantes)
##  [1]  102  137  145  307  731  965 1085 1203 1210 1348 1405 1514 1561 1662 1674
## [16] 1680 1694 1695 1824 1958 2057 2071 2076 2149 2168 2193 2281 2488 2524 2637
## [31] 2653 2666 2849 3058 3155 3169 3227 3231 3243 3245 3294 3319 3418 3461 3493
## [46] 3589 3741 3885 4019 4103 4121
BaseTH[observaciones_faltantes, ]
## # A tibble: 51 × 15
##    SEX     AGE SYSBP CURSMOKE CIGPDAY   BMI DIABETES PREVCHD PREVSTRK PREVHYP
##    <fct> <dbl> <dbl> <fct>      <dbl> <dbl> <fct>    <fct>   <fct>    <fct>  
##  1 0        40  100  0              0  NA   0        0       0        0      
##  2 1        43  110. 1             NA  25.5 0        0       0        0      
##  3 1        49  128. 1             NA  28.2 0        0       0        0      
##  4 0        47  195  1             25  NA   1        0       0        1      
##  5 0        45  108. 0              0  NA   0        0       0        0      
##  6 0        62  242  1             NA  39.9 0        1       1        1      
##  7 0        49  120  1             NA  22.3 0        0       0        0      
##  8 0        64  148  1              3  NA   0        0       0        0      
##  9 0        47  126  0              0  NA   0        0       0        0      
## 10 1        42  122. 1             NA  25.5 0        0       0        0      
## # ℹ 41 more rows
## # ℹ 5 more variables: DEATH <fct>, ANGINA <fct>, HOSPMI <fct>, ANYCHD <fct>,
## #   STROKE <fct>
BaseTH <- BaseTH %>% filter(!row_number() %in% observaciones_faltantes)

y finalmente corroboramos haber eliminado todas las observaciones con datos faltantes

which(!complete.cases(BaseTH))
## integer(0)

/////////// CONSIGNA N1 ///////////

Ahora vamos a calcular la tasa de mortalidad, utilizando la variable DEATH Para eso vamos a crear un objeto que contenga el total de observaciones que tengan “1” en DEATH, y vamos a dividirlo por el total de observaciones * 100 Luego le pedimos que lo imprima para ver el resultado

total_muertes <- sum(BaseTH$DEATH == "1")
tasa_mortalidad <- (total_muertes / nrow(BaseTH)) * 100

print (tasa_mortalidad)
## [1] 34.83915

/////////// CONSIGNA N2 ///////////

Para armar la tabla 1, primero tenemos que definir cuál es la mejor medida para representar a nuestras variables. Sabemos que a las categóricas las representaremos con porcentajes/proporciones, pero para las numéricas debemos definir si representarlas mediante la media y DS; o mediana y RIC según su distribución

Entonces, para las variables numéricas (AGE, SYSBP, CIGPDAY, BMI) lo primero que vamos a hacer, es evaluar variable por variable la distribución de las mismas (normal / no normal) en cada grupo (DEATH = 0/DEATH = 1).

Para esto, vamos a realizar una inspección estadística como también visualización gráfica: 1- En primer lugar, mediante estadísticas descriptivas (media, mediana, desviación estándar, skewness y kurtosis) 2- Ademas, vamos a realizar la prueba de Shapiro-Wilk (H0 distribución normal de los datos) 3- Por último, realizamos histogramas y Q-Q plots, para inspeccionarlas visualmente.

******** variable AGE

## estadística descriptiva
describeBy(BaseTH$AGE, group=BaseTH$DEATH, mat = TRUE) 
##     item group1 vars    n    mean       sd median  trimmed    mad min max range
## X11    1      0    1 2856 47.2437 7.672782     46 46.80184 8.8956  32  69    37
## X12    2      1    1 1527 54.9201 8.218201     56 55.34587 8.8956  34  70    36
##           skew   kurtosis        se
## X11  0.4261453 -0.7169087 0.1435733
## X12 -0.3926248 -0.8280865 0.2103087
## Test de Shapiro-Wilk 
shapiro.test(BaseTH$AGE[BaseTH$DEATH==1])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$AGE[BaseTH$DEATH == 1]
## W = 0.96155, p-value < 2.2e-16
shapiro.test(BaseTH$AGE[BaseTH$DEATH==0])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$AGE[BaseTH$DEATH == 0]
## W = 0.96434, p-value < 2.2e-16
## análisis visual de los resultados con histograma
hist(BaseTH$AGE[BaseTH$DEATH==1]) 

hist(BaseTH$AGE[BaseTH$DEATH==0])

## análisis visual mediante qqplot
qqnorm(BaseTH$AGE[BaseTH$DEATH==1]) 
qqline(BaseTH$AGE[BaseTH$DEATH==1])

qqnorm(BaseTH$AGE[BaseTH$DEATH==0]) 
qqline(BaseTH$AGE[BaseTH$DEATH==0])

AGE –> vamos a considerarla NORMAL (a pesar de test Shapiro con p < 0.05, ya que al ser sensible a tamaños muestrales grandes con valores extremos, puede estar alterado, pero tanto la descripción de los datos como la visualización de los mismos impresiona distribución normal)

******** variable SYSBP

## estadística descriptiva
describeBy(BaseTH$SYSBP, group=BaseTH$DEATH, mat = TRUE)
##     item group1 vars    n     mean       sd median  trimmed     mad  min max
## X11    1      0    1 2856 127.9470 18.56082    125 126.2843 16.3086 85.0 243
## X12    2      1    1 1527 142.0907 25.58694    138 139.9558 23.7216 83.5 295
##     range      skew kurtosis        se
## X11 158.0 1.0175689 1.750167 0.3473106
## X12 211.5 0.8676435 1.095733 0.6547852
## Test de Shapiro-Wilk 
shapiro.test(BaseTH$SYSBP[BaseTH$DEATH==1])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$SYSBP[BaseTH$DEATH == 1]
## W = 0.95689, p-value < 2.2e-16
shapiro.test(BaseTH$SYSBP[BaseTH$DEATH==0])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$SYSBP[BaseTH$DEATH == 0]
## W = 0.94874, p-value < 2.2e-16
## análisis visual de los resultados con histograma
hist(BaseTH$SYSBP[BaseTH$DEATH==1]) 

hist(BaseTH$SYSBP[BaseTH$DEATH==0])

## análisis visual mediante qqplot
qqnorm(BaseTH$SYSBP[BaseTH$DEATH==1])
qqline(BaseTH$SYSBP[BaseTH$DEATH==1])

qqnorm(BaseTH$SYSBP[BaseTH$DEATH==0])
qqline(BaseTH$SYSBP[BaseTH$DEATH==0])

SYSBP –> vamos a considerarla NORMAL (idem anterior, a pesar de test Shapiro con p < 0.05, ya que al ser sensible a tamaños muestrales grandes con valores extremos, puede estar alterado, pero tanto la descripción de los datos como la visualización de los mismos impresiona distribución normal)

******** variable CIGPDAY

## estadística descriptiva
describeBy(BaseTH$CIGPDAY, group=BaseTH$DEATH, mat = TRUE)
##     item group1 vars    n    mean       sd median  trimmed    mad min max range
## X11    1      0    1 2856 8.44958 11.42910      0 6.428259 0.0000   0  70    70
## X12    2      1    1 1527 9.97315 12.78668      1 7.828291 1.4826   0  60    60
##         skew  kurtosis        se
## X11 1.303663 1.2456922 0.2138617
## X12 1.158176 0.6906599 0.3272189
## Test de Shapiro-Wilk 
shapiro.test(BaseTH$CIGPDAY[BaseTH$DEATH==1])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$CIGPDAY[BaseTH$DEATH == 1]
## W = 0.77512, p-value < 2.2e-16
shapiro.test(BaseTH$CIGPDAY[BaseTH$DEATH==0])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$CIGPDAY[BaseTH$DEATH == 0]
## W = 0.75396, p-value < 2.2e-16
## análisis visual de los resultados con histograma
hist(BaseTH$CIGPDAY[BaseTH$DEATH==1])

hist(BaseTH$CIGPDAY[BaseTH$DEATH==0])

## análisis visual mediante qqplot
qqnorm(BaseTH$CIGPDAY[BaseTH$DEATH==1])
qqline(BaseTH$CIGPDAY[BaseTH$DEATH==1])

qqnorm(BaseTH$CIGPDAY[BaseTH$DEATH==0])
qqline(BaseTH$CIGPDAY[BaseTH$DEATH==0])

CIGPDAY –> vamos a considerarla NO NORMAL, ya que tanto el análisis descriptivo como el visual de la misma tiene distribución NO normal.

******** variable BMI

## estadística descriptiva
describeBy(BaseTH$BMI, group=BaseTH$DEATH, mat = TRUE) 
##     item group1 vars    n     mean       sd median  trimmed      mad   min
## X11    1      0    1 2856 25.54720 3.920234  25.11 25.27871 3.528588 16.59
## X12    2      1    1 1527 26.39627 4.365053  26.05 26.14173 3.825108 15.54
##       max range      skew kurtosis         se
## X11 56.80 40.21 1.0390509 3.080554 0.07335553
## X12 51.28 35.74 0.8306954 1.887240 0.11170433
## Test de Shapiro-Wilk 
shapiro.test(BaseTH$BMI[BaseTH$DEATH==1])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$BMI[BaseTH$DEATH == 1]
## W = 0.96537, p-value < 2.2e-16
shapiro.test(BaseTH$BMI[BaseTH$DEATH==0])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$BMI[BaseTH$DEATH == 0]
## W = 0.95423, p-value < 2.2e-16
## análisis visual de los resultados con histograma
hist(BaseTH$BMI[BaseTH$DEATH==1]) 

hist(BaseTH$BMI[BaseTH$DEATH==0])

## análisis visual mediante qqplot
qqnorm(BaseTH$BMI[BaseTH$DEATH==1]) 
qqline(BaseTH$BMI[BaseTH$DEATH==1])

qqnorm(BaseTH$BMI[BaseTH$DEATH==0]) 
qqline(BaseTH$BMI[BaseTH$DEATH==0])

BMI –> vamos a considerarla NORMAL (a pesar de test Shapiro con p < 0.05, ya que al ser sensible a tamaños muestrales grandes con valores extremos, puede estar alterado, pero tanto la descripción de los datos como la visualización de los mismos impresiona distribución normal)

Ahora, conociendo la distribución de las variables, podemos hacer la tabla. Para esto, primero vamos a ponerle labels a las variables y despues armamos la tabla

label(BaseTH$SEX) <- "Sexo masculino, n(%)"
label(BaseTH$AGE) <- "Edad, media(DE)"
label(BaseTH$SYSBP) <- "Presión sistólica, media(DE)"
label(BaseTH$CURSMOKE) <- "Tabaquismo activo, n(%)"
label(BaseTH$CIGPDAY) <- "Cigarrillos por día, mediana(RIC)"
label(BaseTH$BMI) <- "Índice de masa corporal, media(DE)"
label(BaseTH$DIABETES) <- "Antecedente de DBT, n(%)"
label(BaseTH$PREVCHD) <- "Antecedente de Enf Cardíaca, n(%)"
label(BaseTH$PREVSTRK) <- "Antecedente de ACV, n(%)"
label(BaseTH$PREVHYP) <- "Antecedente de HTA, n(%)"

descrTable( DEATH ~ SEX + AGE + SYSBP + CURSMOKE + CIGPDAY + BMI + DIABETES + PREVCHD + PREVSTRK + PREVHYP,
            method=c(SEX=3, AGE=1, SYSBP=1, CURSMOKE=3, CIGPDAY=2, BMI=1, DIABETES=3, PREVCHD=3, PREVSTRK=3, PREVHYP=3),
            BaseTH, hide.no = "0", show.all = TRUE, alpha=0.05,)
## 
## --------Summary descriptives table by 'DEATH'---------
## 
## _______________________________________________________________________________________________ 
##                                         [ALL]              0                1         p.overall 
##                                         N=4383           N=2856           N=1527                
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## Sexo masculino, n(%)                 1923 (43.9%)     1089 (38.1%)     834 (54.6%)     <0.001   
## Edad, media(DE)                      49.9 (8.68)      47.2 (7.67)      54.9 (8.22)     <0.001   
## Presión sistólica, media(DE)          133 (22.3)       128 (18.6)       142 (25.6)     <0.001   
## Tabaquismo activo, n(%)              2142 (48.9%)     1370 (48.0%)     772 (50.6%)      0.109   
## Cigarrillos por día, mediana(RIC)  0.00 [0.00;20.0] 0.00 [0.00;20.0] 1.00 [0.00;20.0]   0.003   
## Índice de masa corporal, media(DE)   25.8 (4.10)      25.5 (3.92)      26.4 (4.37)     <0.001   
## Antecedente de DBT, n(%)             119 (2.72%)       27 (0.95%)       92 (6.02%)     <0.001   
## Antecedente de Enf Cardíaca, n(%)    191 (4.36%)       49 (1.72%)      142 (9.30%)     <0.001   
## Antecedente de ACV, n(%)              29 (0.66%)       7 (0.25%)        22 (1.44%)     <0.001   
## Antecedente de HTA, n(%)             1414 (32.3%)     669 (23.4%)      745 (48.8%)     <0.001   
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

/////////// CONSIGNA N3 ///////////

Para analizar la relación entre la variable SYSBP (variable continua con distribución normal) y DEATH, realizamos un T-test (donde la H0 es que no hay diferencia entre grupos)

t.test(BaseTH$SYSBP~BaseTH$DEATH)
## 
##  Welch Two Sample t-test
## 
## data:  BaseTH$SYSBP by BaseTH$DEATH
## t = -19.082, df = 2403.8, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
## 95 percent confidence interval:
##  -15.59719 -12.69030
## sample estimates:
## mean in group 0 mean in group 1 
##        127.9470        142.0907

Obtenemos una p valor <0.05, por lo cual rechazamos la H0 de no-diferencia entre los grupos, y concluimos que existe diferencia estadísticamente significativa entre los valores de presión arterial sistólica en los fallecidos y no fallecidos. La diferencia entre medias es de 14, con un intervalo de confianza del 95% de entre -15.6 a -12.7

/////////// CONSIGNA N4 ///////////

Ahora vamos a crear una nueva variable llamada CATHYP que será de tipo factor, donde la presión arterial sistólica (variable SYSBP) quede dicotomizada en hipertensos para los valores mayores o iguales a 140, y normotensos para los menores de 140.

BaseTH$CATHYP<-factor(ifelse(BaseTH$SYSBP >=140,"Hipertensos","Normotensos"))

Vamos a analizar si existe relación estadísticamente significativa entre la nueva variable (CATHYP) y DEATH.

Nuestra H0 es que NO existe asociación entre las variables, y la H1 es que existe asociación entre las variables.

Al tratarse de dos variables categóricas, podemos usar test exacto de Fisher o chi2, pero para éste último se debe cumplir el supuesto de N observaciones suficientes. Hacemos una tabla cruzada con ambos métodos y con los valores esperados para decidir cual utilizar.

pas_y_death<- crosstab(BaseTH$CATHYP,BaseTH$DEATH, prop.r = TRUE, prop.c = TRUE, fisher = TRUE, chisq=TRUE, expected =TRUE)

print(pas_y_death)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |         Expected Values | 
## |             Row Percent | 
## |          Column Percent | 
## |-------------------------|
## 
## ========================================
##                  BaseTH$DEATH
## BaseTH$CATHYP         0        1   Total
## ----------------------------------------
## Hipertensos        647      718    1365 
##                  889.4    475.6         
##                   47.4%    52.6%   31.1%
##                   22.7%    47.0%        
## ----------------------------------------
## Normotensos       2209      809    3018 
##                 1966.6   1051.4         
##                   73.2%    26.8%   68.9%
##                   77.3%    53.0%        
## ----------------------------------------
## Total             2856     1527    4383 
##                   65.2%    34.8%        
## ========================================
## 
## Statistics for All Table Factors
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 = 275.4824      d.f. = 1      p <2e-16 
## 
## Pearson's Chi-squared test with Yates' continuity correction 
## ------------------------------------------------------------
## Chi^2 = 274.3473      d.f. = 1      p <2e-16 
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio: 0.3301101 
## 
## Alternative hypothesis: true odds ratio is not equal to 1 
## p <2e-16 
## 95% confidence interval: 0.2881679 0.3779762 
## 
## Alternative hypothesis: true odds ratio is less than 1 
## p <2e-16 
## 95%s confidence interval: % 0 0.369985 
## 
## Alternative hypothesis: true odds ratio is greater than 1 
## p = 1 
## 95%s confidence interval: % 0.2944444 Inf 
## 
##         Minimum expected frequency: 475.5544

Como la frecuencia mínima esperada es de 475 (en todas las celdas hay N suficiente), se cumple el supuesto de suficientes observaciones y podemos utilizar cualquier método.

Ambos test de hipótesis arrojan un p valor <0.05, por lo cual podemos rechazar la H0 de NO asociación entre las variables y concluir que existe asociación estadísticamente significativa entre pacientes hipertensos y el fallecimiento.

Ahora, vamos a ver la magnitud de la asociación entre estas variables. Para eso, vamos a calcular el RR de las mismas (ya que se trata de un estudio de cohorte).

epi.2by2(table(BaseTH$CATHYP,BaseTH$DEATH)[c("Hipertensos","Normotensos"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +          718          647       1365     52.60 (49.91 to 55.28)
## Exposed -          809         2209       3018     26.81 (25.23 to 28.42)
## Total             1527         2856       4383     34.84 (33.43 to 36.27)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 1.96 (1.82, 2.12)
## Inc odds ratio                                 3.03 (2.65, 3.46)
## Attrib risk in the exposed *                   25.79 (22.71, 28.88)
## Attrib fraction in the exposed (%)            49.04 (44.93, 52.84)
## Attrib risk in the population *                8.03 (5.92, 10.15)
## Attrib fraction in the population (%)         23.06 (20.09, 25.92)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 275.482 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

RR: 1.96 (IC95% 1.82-2.12). Esto se puede interpretar como: el riesgo de muerte en los pacientes hipertensos es 96% mayor que en el grupo de pacientes normotensos, y se puede afirmar con un 95% de certeza que este aumento del riesgo se encuentra entre 82 y 112%.

/////////// CONSIGNA N5 ///////////

Para evaluar si la variable CURSMOKE se comporta como potencial confundidor entre la variable CATHYP (exposición) y DEATH (outcome), primero vamos a evaluar si existe asociación con cada una por separado.

Al tratarse de todas variables categóricas, vamos a usar un test de Fisher:

## asociación entre CATHYP y CURSMOKE. H0 las variables NO están asociadas
pas_y_tbq <- crosstab(BaseTH$CATHYP,BaseTH$CURSMOKE, prop.r = TRUE, prop.c = TRUE, fisher = TRUE, expected =TRUE)

print(pas_y_tbq)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |         Expected Values | 
## |             Row Percent | 
## |          Column Percent | 
## |-------------------------|
## 
## ========================================
##                  Tabaquismo activo, n(%)
## BaseTH$CATHYP         0        1   Total
## ----------------------------------------
## Hipertensos        820      545    1365 
##                  697.9    667.1         
##                   60.1%    39.9%   31.1%
##                   36.6%    25.4%        
## ----------------------------------------
## Normotensos       1421     1597    3018 
##                 1543.1   1474.9         
##                   47.1%    52.9%   68.9%
##                   63.4%    74.6%        
## ----------------------------------------
## Total             2241     2142    4383 
##                   51.1%    48.9%        
## ========================================
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio: 1.690699 
## 
## Alternative hypothesis: true odds ratio is not equal to 1 
## p = 1.52e-15 
## 95% confidence interval: 1.482083 1.929632 
## 
## Alternative hypothesis: true odds ratio is less than 1 
## p = 1 
## 95%s confidence interval: % 0 1.889629 
## 
## Alternative hypothesis: true odds ratio is greater than 1 
## p = 8.93e-16 
## 95%s confidence interval: % 1.513201 Inf 
## 
##         Minimum expected frequency: 667.0842
## asociación entre CURSMOKE y DEATH. H0 las variables NO están asociadas
tbq_y_death <- crosstab(BaseTH$DEATH,BaseTH$CURSMOKE, prop.r = TRUE, prop.c = TRUE, fisher = TRUE, expected =TRUE)

print(tbq_y_death)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |         Expected Values | 
## |             Row Percent | 
## |          Column Percent | 
## |-------------------------|
## 
## =======================================
##                 Tabaquismo activo, n(%)
## BaseTH$DEATH         0        1   Total
## ---------------------------------------
## 0                1486     1370    2856 
##                1460.3   1395.7         
##                  52.0%    48.0%   65.2%
##                  66.3%    64.0%        
## ---------------------------------------
## 1                 755      772    1527 
##                 780.7    746.3         
##                  49.4%    50.6%   34.8%
##                  33.7%    36.0%        
## ---------------------------------------
## Total            2241     2142    4383 
##                  51.1%    48.9%        
## =======================================
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio: 1.109067 
## 
## Alternative hypothesis: true odds ratio is not equal to 1 
## p = 0.106 
## 95% confidence interval: 0.9775007 1.258392 
## 
## Alternative hypothesis: true odds ratio is less than 1 
## p = 0.952 
## 95%s confidence interval: % 0 1.233463 
## 
## Alternative hypothesis: true odds ratio is greater than 1 
## p = 0.0547 
## 95%s confidence interval: % 0.9972193 Inf 
## 
##         Minimum expected frequency: 746.2546

En primera instancia, al no encontrar asociación estadística entre CURSMOKE y DEATH, no sería un potencial confundidor.

Para ver si es modificador de efecto, vamos a correr un test Mantel Hanzel estratificando por CURSMOKE –> Vamos a testear la hipótesis H0 de homogeneidad entre los estratos.

epi.2by2(table(BaseTH$CATHYP,BaseTH$DEATH,BaseTH$CURSMOKE)[c("Hipertensos","Normotensos"),c("1","0"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +          718          647       1365     52.60 (49.91 to 55.28)
## Exposed -          809         2209       3018     26.81 (25.23 to 28.42)
## Total             1527         2856       4383     34.84 (33.43 to 36.27)
## 
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio (crude)                         1.96 (1.82, 2.12)
## Inc risk ratio (M-H)                           2.01 (1.85, 2.17)
## Inc risk ratio (crude:M-H)                     0.98
## Inc odds ratio (crude)                         3.03 (2.65, 3.46)
## Inc odds ratio (M-H)                           3.14 (2.75, 3.60)
## Inc odds ratio (crude:M-H)                     0.96
## Attrib risk in the exposed (crude) *           25.79 (22.71, 28.88)
## Attrib risk in the exposed (M-H) *             26.48 (21.36, 31.61)
## Attrib risk (crude:M-H)                        0.97
## -------------------------------------------------------------------
##  M-H test of homogeneity of IRRs: chi2(1) = 2.311 Pr>chi2 = 0.128
##  M-H test of homogeneity of ORs: chi2(1) = 0.500 Pr>chi2 = 0.480
##  Test that M-H adjusted OR = 1:  chi2(1) = 286.680 Pr>chi2 = <0.001
##  Wald confidence limits
##  M-H: Mantel-Haenszel; CI: confidence interval
##  * Outcomes per 100 population units

Vemos que el test de MH es no significativo (p>0.05) por lo tanto no rechaza la H0 de homogeneidad, de todas maneras vamos a realizar un análisis estratificado para corroborar estos resultados.

Para el estrato de CURSMOKE si (1)

epi.2by2(table(BaseTH$CATHYP[BaseTH$CURSMOKE==1],BaseTH$DEATH[BaseTH$CURSMOKE==1])[c("Hipertensos","Normotensos"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +          302          243        545     55.41 (51.13 to 59.64)
## Exposed -          470         1127       1597     29.43 (27.20 to 31.73)
## Total              772         1370       2142     36.04 (34.00 to 38.12)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 1.88 (1.69, 2.10)
## Inc odds ratio                                 2.98 (2.44, 3.64)
## Attrib risk in the exposed *                   25.98 (21.25, 30.72)
## Attrib fraction in the exposed (%)            46.89 (40.89, 52.28)
## Attrib risk in the population *                6.61 (3.59, 9.63)
## Attrib fraction in the population (%)         18.34 (14.68, 21.85)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 119.001 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Para el estrato de CURSMOKE no (0)

epi.2by2(table(BaseTH$CATHYP[BaseTH$CURSMOKE==0],BaseTH$DEATH[BaseTH$CURSMOKE==0])[c("Hipertensos","Normotensos"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +          416          404        820     50.73 (47.25 to 54.21)
## Exposed -          339         1082       1421     23.86 (21.66 to 26.16)
## Total              755         1486       2241     33.69 (31.73 to 35.69)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 2.13 (1.90, 2.39)
## Inc odds ratio                                 3.29 (2.74, 3.95)
## Attrib risk in the exposed *                   26.88 (22.80, 30.95)
## Attrib fraction in the exposed (%)            52.98 (47.26, 58.08)
## Attrib risk in the population *                9.83 (6.88, 12.79)
## Attrib fraction in the population (%)         29.19 (24.42, 33.66)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 168.108 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Cómo vemos, si bien los RR de muerte en función de la categoría de HTA no son iguales, estos varían muy ligeramente (RR para el tbq actual 1.9 y para el no tbq 2.1) y los intervalos de confianza se superponen ampliamente.

/////////// CONSIGNA N6 ///////////

Para evaluar si existe asociación entre la variable DIABETES y la variable DEATH, nuevamente como ambas son variables categóricas, podemos usar test exacto de Fisher o chi2, pero para éste último se debe cumplir el supuesto de N observaciones suficientes.

Nuestra H0 es que no existe asociación entre ambas variables. Hacemos una tabla cruzada con ambos métodos para decidir cual utilizar.

diabetes_y_death <- crosstab(BaseTH$DIABETES,BaseTH$DEATH, prop.r = TRUE, prop.c = TRUE, chisq=TRUE, fisher = TRUE, expected =TRUE)

print (diabetes_y_death)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |         Expected Values | 
## |             Row Percent | 
## |          Column Percent | 
## |-------------------------|
## 
## ===================================================
##                             BaseTH$DEATH
## Antecedente de DBT, n(%)         0        1   Total
## ---------------------------------------------------
## 0                            2829     1435    4264 
##                            2778.5   1485.5         
##                              66.3%    33.7%   97.3%
##                              99.1%    94.0%        
## ---------------------------------------------------
## 1                              27       92     119 
##                              77.5     41.5         
##                              22.7%    77.3%    2.7%
##                               0.9%     6.0%        
## ---------------------------------------------------
## Total                        2856     1527    4383 
##                              65.2%    34.8%        
## ===================================================
## 
## Statistics for All Table Factors
## 
## Pearson's Chi-squared test 
## ------------------------------------------------------------
## Chi^2 = 97.19585      d.f. = 1      p <2e-16 
## 
## Pearson's Chi-squared test with Yates' continuity correction 
## ------------------------------------------------------------
## Chi^2 = 95.28227      d.f. = 1      p <2e-16 
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio: 6.714384 
## 
## Alternative hypothesis: true odds ratio is not equal to 1 
## p <2e-16 
## 95% confidence interval: 4.308412 10.77926 
## 
## Alternative hypothesis: true odds ratio is less than 1 
## p = 1 
## 95%s confidence interval: % 0 9.995072 
## 
## Alternative hypothesis: true odds ratio is greater than 1 
## p <2e-16 
## 95%s confidence interval: % 4.597466 Inf 
## 
##         Minimum expected frequency: 41.45859

Vemos que podemos utilizar cualquiera de los dos métodos, por lo que corremos la siguiente

Al tener un valor de p <0.05, descartamos la H0 de NO asociación, y por ende concluimos que EXISTE ASOCIACIÓN estadísticamente significativa entre las variables.

Ahora para CUANTIFICAR esa asociación, calculamos un RR

epi.2by2(table(BaseTH$DIABETES,BaseTH$DEATH)[c("1","0"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +           92           27        119     77.31 (68.73 to 84.48)
## Exposed -         1435         2829       4264     33.65 (32.24 to 35.09)
## Total             1527         2856       4383     34.84 (33.43 to 36.27)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 2.30 (2.07, 2.55)
## Inc odds ratio                                 6.72 (4.35, 10.36)
## Attrib risk in the exposed *                   43.66 (36.00, 51.31)
## Attrib fraction in the exposed (%)            56.47 (51.60, 60.85)
## Attrib risk in the population *                1.19 (-0.81, 3.19)
## Attrib fraction in the population (%)         3.40 (2.55, 4.25)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 97.196 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

El riesgo relativo es de 2.3, con un intervalo de confianza del 95% de 2.07 a 2.55. Lo cual se puede interpretar como: El riesgo de fallecimiento en los pacientes diabéticos es 130% mayor que en el grupo de pacientes no diabéticos, y se puede afirmar con un 95% de certeza que este aumento del riesgo se encuentra entre 107 y 155%.

/////////// CONSIGNA N7 ///////////

Para evaluar si la variable SEX se comporta como potencial confundidor entre la variable DIABETES (exposición) y DEATH (outcome), nuevamente vamos primero a evaluar si existe asociación con cada una por separado.

Para calcular esto, primero analizamos la asociación entre la variable SEXO con la exposición (diabetes) y luego con el outcome (death) mediante un test de Fisher.

## asociación entre SEX y DIABETES, H0 las variables NO están asociadas
sex_y_diabetes <- crosstab(BaseTH$SEX,BaseTH$DIABETES, prop.r = TRUE, prop.c = TRUE, fisher = TRUE, expected = TRUE )

print (sex_y_diabetes)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |         Expected Values | 
## |             Row Percent | 
## |          Column Percent | 
## |-------------------------|
## 
## ==============================================
##                         Antecedente de DBT, n(%)
## Sexo masculino, n(%)         0       1   Total
## ----------------------------------------------
## 0                        2399      61    2460 
##                        2393.2    66.8         
##                          97.5%    2.5%   56.1%
##                          56.3%   51.3%        
## ----------------------------------------------
## 1                        1865      58    1923 
##                        1870.8    52.2         
##                          97.0%    3.0%   43.9%
##                          43.7%   48.7%        
## ----------------------------------------------
## Total                    4264     119    4383 
##                          97.3%    2.7%        
## ==============================================
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio: 1.223003 
## 
## Alternative hypothesis: true odds ratio is not equal to 1 
## p = 0.303 
## 95% confidence interval: 0.8342929 1.791266 
## 
## Alternative hypothesis: true odds ratio is less than 1 
## p = 0.88 
## 95%s confidence interval: % 0 1.689055 
## 
## Alternative hypothesis: true odds ratio is greater than 1 
## p = 0.161 
## 95%s confidence interval: % 0.8850361 Inf 
## 
##         Minimum expected frequency: 52.21013
## asociación entre SEX y DEATH, H0 las variables NO están asociadas
sex_y_death <- crosstab(BaseTH$SEX,BaseTH$DEATH, prop.r = TRUE, prop.c = TRUE, fisher = TRUE, expected = TRUE )

print (sex_y_death)
##    Cell Contents 
## |-------------------------|
## |                   Count | 
## |         Expected Values | 
## |             Row Percent | 
## |          Column Percent | 
## |-------------------------|
## 
## =============================================
##                         BaseTH$DEATH
## Sexo masculino, n(%)        0       1   Total
## ---------------------------------------------
## 0                       1767     693    2460 
##                         1603     857         
##                         71.8%   28.2%   56.1%
##                         61.9%   45.4%        
## ---------------------------------------------
## 1                       1089     834    1923 
##                         1253     670         
##                         56.6%   43.4%   43.9%
##                         38.1%   54.6%        
## ---------------------------------------------
## Total                   2856    1527    4383 
##                         65.2%   34.8%        
## =============================================
## 
##  
## Fisher's Exact Test for Count Data
## ------------------------------------------------------------
## Sample estimate odds ratio: 1.952441 
## 
## Alternative hypothesis: true odds ratio is not equal to 1 
## p <2e-16 
## 95% confidence interval: 1.718192 2.219195 
## 
## Alternative hypothesis: true odds ratio is less than 1 
## p = 1 
## 95%s confidence interval: % 0 2.17464 
## 
## Alternative hypothesis: true odds ratio is greater than 1 
## p <2e-16 
## 95%s confidence interval: % 1.753253 Inf 
## 
##         Minimum expected frequency: 669.9569

Vemos que tiene asociación estadísticamente significativa con la muerte, pero no con la diabetes, así que no parece ser un confundidor.

Para evaluar si la misma de todas maneras puede actuar como modificadora de efecto, vamos a correr un test Mantel Haenzel ajustando por SEX –> Vamos a testear la hipótesis H0 de homogeneidad entre los estratos.

epi.2by2(table(BaseTH$DIABETES,BaseTH$DEATH,BaseTH$SEX)[c("1","0"),c("1","0"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +           92           27        119     77.31 (68.73 to 84.48)
## Exposed -         1435         2829       4264     33.65 (32.24 to 35.09)
## Total             1527         2856       4383     34.84 (33.43 to 36.27)
## 
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio (crude)                         2.30 (2.07, 2.55)
## Inc risk ratio (M-H)                           2.25 (2.03, 2.49)
## Inc risk ratio (crude:M-H)                     1.02
## Inc odds ratio (crude)                         6.72 (4.35, 10.36)
## Inc odds ratio (M-H)                           6.70 (4.31, 10.40)
## Inc odds ratio (crude:M-H)                     1.00
## Attrib risk in the exposed (crude) *           43.66 (36.00, 51.31)
## Attrib risk in the exposed (M-H) *             42.91 (-34.42, 120.23)
## Attrib risk (crude:M-H)                        1.02
## -------------------------------------------------------------------
##  M-H test of homogeneity of IRRs: chi2(1) = 10.580 Pr>chi2 = 0.001
##  M-H test of homogeneity of ORs: chi2(1) = 0.388 Pr>chi2 = 0.533
##  Test that M-H adjusted OR = 1:  chi2(1) = 95.352 Pr>chi2 = <0.001
##  Wald confidence limits
##  M-H: Mantel-Haenszel; CI: confidence interval
##  * Outcomes per 100 population units

En este test rechazamos la H0 de homogeneidad entre estratos. Por lo cual vamos a realizar una estratificación para confirmar estos resultados.

epi.2by2(table(BaseTH$DIABETES[BaseTH$SEX==1],BaseTH$DEATH[BaseTH$SEX==1])[c("1","0"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +           47           11         58     81.03 (68.59 to 90.13)
## Exposed -          787         1078       1865     42.20 (39.94 to 44.48)
## Total              834         1089       1923     43.37 (41.14 to 45.62)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 1.92 (1.68, 2.20)
## Inc odds ratio                                 5.85 (3.02, 11.36)
## Attrib risk in the exposed *                   38.84 (28.50, 49.17)
## Attrib fraction in the exposed (%)            47.93 (40.38, 54.52)
## Attrib risk in the population *                1.17 (-1.98, 4.32)
## Attrib fraction in the population (%)         2.70 (1.70, 3.69)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 34.543 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units
epi.2by2(table(BaseTH$DIABETES[BaseTH$SEX==0],BaseTH$DEATH[BaseTH$SEX==0])[c("1","0"),c("1","0")], method = "cohort.count", conf.level = 0.95, outcome = "as.columns")
##              Outcome +    Outcome -      Total                 Inc risk *
## Exposed +           45           16         61     73.77 (60.93 to 84.20)
## Exposed -          648         1751       2399     27.01 (25.24 to 28.84)
## Total              693         1767       2460     28.17 (26.40 to 29.99)
## 
## Point estimates and 95% CIs:
## -------------------------------------------------------------------
## Inc risk ratio                                 2.73 (2.32, 3.22)
## Inc odds ratio                                 7.60 (4.27, 13.54)
## Attrib risk in the exposed *                   46.76 (35.58, 57.94)
## Attrib fraction in the exposed (%)            63.38 (56.88, 68.91)
## Attrib risk in the population *                1.16 (-1.35, 3.67)
## Attrib fraction in the population (%)         4.12 (2.69, 5.52)
## -------------------------------------------------------------------
## Uncorrected chi2 test that OR = 1: chi2(1) = 64.278 Pr>chi2 = <0.001
## Fisher exact test that OR = 1: Pr>chi2 = <0.001
##  Wald confidence limits
##  CI: confidence interval
##  * Outcomes per 100 population units

Y acá efectivamente vemos que en el grupo de hombres, el IRR es de 1.93 mientras que en mujeres de 2.73. Concluimos asi que no es una confundidor pero si un modificador de efecto.

/////////// CONSIGNA N8 ///////////

Primero genero las categorías de la variable, luego la transformo en factor y reviso como quedaron los niveles de las categorias.

BaseTH <- BaseTH %>% mutate(OBESE = case_when( BMI <= 25 ~ "normopeso", BMI > 25 & BMI <= 30 ~ "sobrepeso",BMI > 30 ~ "obesidad" ))

BaseTH$OBESE<-as.factor(BaseTH$OBESE)  
is.factor(BaseTH$OBESE)
## [1] TRUE
levels(BaseTH$OBESE)
## [1] "normopeso" "obesidad"  "sobrepeso"

/////////// CONSIGNA N9 ///////////

Para comparar la presion sistólica (variable continua) en las 3 categorías de OBESE podemos utilizar ANOVA o su alternativa no parámetrica Kruskal-Wallis. Para usar ANOVA, se deben cumplir los siguientes 3 supuestos:

describeBy(BaseTH$SYSBP, group=BaseTH$OBESE, mat = TRUE) 
##     item    group1 vars    n     mean       sd median  trimmed     mad  min max
## X11    1 normopeso    1 1979 126.5915 20.35418    123 124.2798 16.3086 83.5 244
## X12    2  obesidad    1  572 144.9327 24.28619    141 142.9607 20.7564 93.0 295
## X13    3 sobrepeso    1 1832 135.8968 21.51270    132 133.7664 19.2738 83.5 243
##     range      skew kurtosis        se
## X11 160.5 1.3917252 3.131235 0.4575418
## X12 202.0 1.0502805 2.590104 1.0154566
## X13 159.5 0.9813999 1.116899 0.5026112
shapiro.test(BaseTH$SYSBP[BaseTH$OBESE == "normopeso"])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$SYSBP[BaseTH$OBESE == "normopeso"]
## W = 0.91221, p-value < 2.2e-16
shapiro.test(BaseTH$SYSBP[BaseTH$OBESE == "obesidad"])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$SYSBP[BaseTH$OBESE == "obesidad"]
## W = 0.94804, p-value = 2.748e-13
shapiro.test(BaseTH$SYSBP[BaseTH$OBESE == "sobrepeso"])
## 
##  Shapiro-Wilk normality test
## 
## data:  BaseTH$SYSBP[BaseTH$OBESE == "sobrepeso"]
## W = 0.94499, p-value < 2.2e-16
hist(BaseTH$SYSBP[BaseTH$OBESE == "normopeso"])

hist(BaseTH$SYSBP[BaseTH$OBESE == "obesidad"])

hist(BaseTH$SYSBP[BaseTH$OBESE == "sobrepeso"])

qqnorm(BaseTH$SYSBP[BaseTH$OBESE == "normopeso"])
qqline(BaseTH$SYSBP[BaseTH$OBESE == "normopeso"])

qqnorm(BaseTH$SYSBP[BaseTH$OBESE == "obesidad"])
qqline(BaseTH$SYSBP[BaseTH$OBESE == "obesidad"])

qqnorm(BaseTH$SYSBP[BaseTH$OBESE == "sobrepeso"])
qqline(BaseTH$SYSBP[BaseTH$OBESE == "sobrepeso"])

bartlett.test(BaseTH$SYSBP,BaseTH$OBESE)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  BaseTH$SYSBP and BaseTH$OBESE
## Bartlett's K-squared = 29.618, df = 2, p-value = 3.703e-07
leveneTest(BaseTH$SYSBP,BaseTH$OBESE)
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value    Pr(>F)    
## group    2  16.376 8.211e-08 ***
##       4380                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Habiendo definido que la distribución es de tipo NO normal y de que no existe homogeneidad de varianzas (ambos test rechazan la H0 de homogeneidad), no podemos utilizar ANOVA, por lo que optamos por un test no paramétrico como Kruskal Wallis (con la H0 de que no existe diferencia entre las categorías y Ha de que al menos uno de los grupos es diferente)

kruskal.test(SYSBP~OBESE, data = BaseTH)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  SYSBP by OBESE
## Kruskal-Wallis chi-squared = 412.98, df = 2, p-value < 2.2e-16

Obtenemos una p <0.05, por lo que rechazamos la H0 de que no existe diferencia entre las categorías.

Sabemos que al menos una de las categorias es diferentes, y para eso correremos un wilcoxon-test de a pares (para identificar cual o cuales categorías de obesidad tienen diferencia estadísticamente significativa), hacemos una corrección de bonferroni ya que estamos realizando comparaciones múltiples.

La H0 del wilcoxon-test es que ambos grupos son iguales.

pairwise.wilcox.test(BaseTH$SYSBP,BaseTH$OBESE,"bonferroni")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  BaseTH$SYSBP and BaseTH$OBESE 
## 
##           normopeso obesidad
## obesidad  <2e-16    -       
## sobrepeso <2e-16    <2e-16  
## 
## P value adjustment method: bonferroni

Obtenemos p valor <0.05 en todas las comparaciones, por lo que descartamos la H0 de igualdad entre grupos y cocluimos que existe diferencia estadísticamente significativa en los valores de presión arterial sistólicas para todas las categorías de obesidad.

/////////// CONSIGNA N10 /////////// –> tabla realizada en el word

/////////// CONSIGNA N11 ///////////

+++++ SEX según DEATH

# Primero le otorgamos labels a las observaciones de cada variable, para que queden visualmente mas facil de entender en los gráficos
BaseTH <- BaseTH %>%
  mutate(
    SEX = factor(SEX, labels = c("Mujeres", "Hombres")),
    CURSMOKE = factor(CURSMOKE, labels = c("No tabaquista", "Tabaquista")),
    DIABETES = factor(DIABETES, labels = c("No diabéticos", "Diabéticos")),
    CATHYP = factor(CATHYP, labels = c("Hipertensos","Normotensos")),
    DEATH = factor(DEATH, labels = c("No fallecidos", "Fallecidos"))
  )
# Creamos el gráfico con colores específicos
ggplot(data = BaseTH, aes(x = DEATH, fill = SEX)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("pink", "skyblue")) +
  labs(y = "Proporción por cada género", x = "Fallecimientos", fill = "Género", title="Proporciones de sexo femenino y masculino según sobrevida") +
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

+++++ SYSBP según DEATH

ggplot(data = BaseTH, aes(x = DEATH, y = SYSBP)) +
    geom_boxplot(fill = c("grey")) +
    labs(x = "Fallecimientos", y = "Presión sistólica (mmHg)", title = "Distribución de presión sistólica según sobrevida") +
    theme_minimal() +
    theme(plot.title = element_text(hjust = 0.5))

+++++ CURSMOKE según DEATH

# Creamos el gráfico con colores específicos
ggplot(data = BaseTH, aes(x = DEATH, fill = CURSMOKE)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgreen", "red")) +
  labs(y = "Proporción de tabaquistas", x = "Fallecimientos", fill = "Tabaquismo", title= "Proporción de tabaquismo según sobrevida") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

+++++ BMI según DEAT

ggplot(data = BaseTH, aes(x = BMI, fill = factor(DEATH))) +
  geom_density(alpha = 0.5) +
  scale_fill_manual(values = c("blue", "red")) +
  labs(x = "IMB", y = "Fallecimientos", title = "Distribución de IMC según supervivencia") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

+++++ DIABETES según DEATH

# Creamos el gráfico con colores específicos
ggplot(data = BaseTH, aes(x = DEATH, fill = DIABETES)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgreen", "red")) +
  labs(y = "Proporción de diabéticos", x = "Fallecimientos", fill = "Diabetes", title="Proporción de diabetes según sobrevida") +
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))

+++++++ CATHYP según DEATH

# Creamos el gráfico con colores específicos
ggplot(data = BaseTH, aes(x = DEATH, fill = CATHYP)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = c("lightgreen", "red")) +
  labs(y = "Proporción de hipertensos", x = "Fallecimientos", fill = "Hipertensión arterial", title="Proporción de hipertensión arterial según sobrevida") +
  theme_minimal()+
  theme(plot.title = element_text(hjust = 0.5))