En esta sesión veremos diversas herramientas disponibles en R para la realización de contrastes de hipótesis y el cálculo de intervalos de confianza.
Como en anteriores ocasiones, utilizaremos los datos del estudio de Telde:
library(openxlsx)
setwd("c:/Users/aulas/Downloads/")
telde = read.xlsx("endocrino.xlsx")
Para determinar si el sexo de una persona se asocia con el hecho de que esa persona padezca o no HTA comenzamos por la construcción de una tabla cruzada de ambas variables:
telde$SEXO=ordered(telde$SEXO,levels=c(1,0),labels=c("Hombre","Mujer"))
telde$HTA=ordered(telde$HTA_OMS,levels=c(1,0),labels=c("HTA+","HTA-"))
tb=with(telde,table(HTA,SEXO))
tb
## SEXO
## HTA Hombre Mujer
## HTA+ 156 168
## HTA- 292 414
prop.table(tb,2)
## SEXO
## HTA Hombre Mujer
## HTA+ 0.3482143 0.2886598
## HTA- 0.6517857 0.7113402
Esta tabla indica que un 34.8% de los hombres son hipertensos, frente a un 28.9% de las mujeres, lo que podría interpretarse como que el ser hombre es de algún modo un factor de riesgo de hipertensión. Podemos evaluar la asociación entre sexo y HTA mediante la odds-ratio:
tb[1,1]*tb[2,2]/(tb[1,2]*tb[2,1])
## [1] 1.316536
cuyo valor (mayor que la unidad) está efectivamente acorde con la presunción de que ser hombre incrementa el riesgo de HTA.
En cualquier caso la diferencia entre un 34.8% y un 28.9% no parece excesivamente elevada. Dado que los datos con los que estamos trabajando constituyen una muestra elegida al azar entre la población adulta de Telde cabe hacerse la siguiente pregunta: si la muestra hubiese sido otra ¿se habría observado la misma relación entre estos porcentajes (esto es mayor porcentaje de hipertensos entre las mujeres que entre los hombres), o quizás estos porcentajes podrían llegar incluso a invertirse? En otras palabras, la pregunta que nos hacemos es: ¿la diferencia observada puede explicarse simplemente por la variabilidad aleatoria presente en el muestreo, o necesariamente se debe a que en realidad la hipertensión es más común entre los hombres que entre las mujeres de Telde?
Para responder a esta pregunta debe utilizarse el test de la chi-cuadrado:
\[\frac{\pi_M}{\pi_H}=1\]
\[\frac{\pi_M}{\pi_H}\neq 1\]
La realización del test de la chi-cuadrado en R es muy simple: se utiliza la función chisq.test
a la que se le pasa como argumento la tabla con los datos disponibles:
chisq.test(tb)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tb
## X-squared = 3.8924, df = 1, p-value = 0.0485
Si deseamos calcular las frecuencias esperadas bajo la hipótesis nula (igual prevalencia de HTA en hombres y mujeres) podemos utilizar:
chisq.test(tb)$expected
## SEXO
## HTA Hombre Mujer
## HTA+ 140.9243 183.0757
## HTA- 307.0757 398.9243
Para responder a esta pregunta podemos proceder del mismo modo que con el sexo:
telde$INSTRUCCION=ordered(telde$INSTRUCCION,levels=c("Sin estudios","Primer grado",
"Segundo grado","Tercer grado"))
tb=with(telde,table(HTA,INSTRUCCION))
tb
## INSTRUCCION
## HTA Sin estudios Primer grado Segundo grado Tercer grado
## HTA+ 6 167 130 21
## HTA- 7 191 379 129
prop.table(tb,2)
## INSTRUCCION
## HTA Sin estudios Primer grado Segundo grado Tercer grado
## HTA+ 0.4615385 0.4664804 0.2554028 0.1400000
## HTA- 0.5384615 0.5335196 0.7445972 0.8600000
Esta última tabla muestra que a medida que aumenta el nivel de instrucción aparentemente disminuye la probabilidad de padecer hipertensión.
Aplicamos el test de la chi-cuadrado:
chisq.test(tb)
## Warning in chisq.test(tb): Chi-squared approximation may be incorrect
##
## Pearson's Chi-squared test
##
## data: tb
## X-squared = 69.084, df = 3, p-value = 6.705e-15
Observemos que R nos muestra una advertencia relativa a que la aproximación chi-cuadrado puede ser incorrecta. Ello es debido a que las condiciones “técnicas” de realización del test requieren que ninguna frecuencia esperada sea menor que 1 y que al menos el 80 % de las frecuencias esperadas son mayores que 5. Las frecuencias esperadas en este caso son:
chisq.test(tb)$expected
## Warning in chisq.test(tb): Chi-squared approximation may be incorrect
## INSTRUCCION
## HTA Sin estudios Primer grado Segundo grado Tercer grado
## HTA+ 4.08932 112.6136 160.1126 47.18447
## HTA- 8.91068 245.3864 348.8874 102.81553
En este caso hay una única frecuencia esperada menor que 5; en total hay 8 valores esperados de los que 7, esto es el 87.5%, son mayores que 5. Por tanto podemos confiar en el resultado del test y concluir que la asociación entre el nivel de instrucción y la HTA es significativa, esto es, no puede explicarse simplemente por efecto del azar.
En caso de que no se hubiesen cumplido las condiciones para la realización del test, se debería proceder a agrupar categorías; en este caso, como la categoría que puede causar problemas es la categoría “Sin estudios”, ésta podría fusionarse con la categoría “Primer grado”:
telde$INSTRUCCION2=ifelse((telde$INSTRUCCION=="Sin estudios")|(telde$INSTRUCCION=="Primer grado"),"Sin Estudios/Primer grado",as.character(telde$INSTRUCCION))
tb=with(telde,table(HTA,INSTRUCCION2))
tb
## INSTRUCCION2
## HTA Segundo grado Sin Estudios/Primer grado Tercer grado
## HTA+ 130 173 21
## HTA- 379 198 129
prop.table(tb,2)
## INSTRUCCION2
## HTA Segundo grado Sin Estudios/Primer grado Tercer grado
## HTA+ 0.2554028 0.4663073 0.1400000
## HTA- 0.7445972 0.5336927 0.8600000
chisq.test(tb)
##
## Pearson's Chi-squared test
##
## data: tb
## X-squared = 69.082, df = 2, p-value = 9.976e-16
El índice de masa corporal se calcula como el peso (en kg) partido por la talla (en m) elevada al cuadrado. Podemos calcular el IMC de las personas de la muestra de Telde mediante:
telde$IMC=telde$PESO/(telde$TALLA/100)^2
La siguiente tabla nos da el valor medio de IMC según las personas sean o no hipertensas:
aggregate (IMC~HTA, telde, mean)
## HTA IMC
## 1 HTA+ 31.09140
## 2 HTA- 26.92119
Gráficamente:
library(ULPGCmisc)
summarize(telde$IMC,by=telde$HTA, showSummary = FALSE)
Vemos, pues, que el IMC medio es mayor en los hipertensos (31.10) que en los normotensos (26.92). Igual que en el caso anterior: ¿podría esta diferencia ser efecto de la variabilidad aleatoria del muestreo o la diferencia es tan grande que solo puede explicarse porque en realidad un mayor IMC se asocia con riesgo de HTA?
Para responder a esta pregunta podemos utilizar el test de la de Student:
La hipótesis nula es que la hipertensión no se asocia con el IMC; ello significaría que el valor medio de IMC sería el mismo en el grupo de hipertensos y en el grupo de normotensos. Si llamamos \(\mu_H\) al valor medio de IMC en el primer grupo y \(\mu_N\) en el segundo, la hipótesis nula especifica que: \[\mu_H=\mu_N\]
La hipótesis alternativa es que el valor medio de IMC difiere entre hipertensos y normotensos; podemos ser más específicos y considerar que la hipótesis alternativa es: \[\mu_H>\mu_N\]
esto es que el valor medio de IMC es estrictamente mayor en hipertensos que en normotensos. Para resolver este contraste ejecutamos la siguiente función en R:
t.test(IMC~HTA,data=telde,alternative="greater")
##
## Welch Two Sample t-test
##
## data: IMC by HTA
## t = 12.909, df = 565.08, p-value < 2.2e-16
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 3.637958 Inf
## sample estimates:
## mean in group HTA+ mean in group HTA-
## 31.09140 26.92119
El p-valor es muy pequeño (<0.0001) lo que indica que hay evidencia suficiente para asegurar que el valor de IMC en el grupo de hipertensos es mayor que en el grupo de los normotensos; la diferencia observada no puede explicarse por el mero efecto del azar.
La salida de R nos informa además que con una confianza del 95% el valor medio de IMC en hipertensos es al menos 3.64 unidades mayor que en normotensos.
Si queremos obtener un intervalo de confianza al 95% para la diferencia en el valor medio de IMC entre hipertensos y normotensos utilizamos la sintaxis anterior sin especificar alternativa; asimismo podemos indicar fácilmente el nivel de confianza deseado:
t.test(IMC~HTA,data=telde, conf=0.99)$conf.int
## [1] 3.335257 5.005155
## attr(,"conf.level")
## [1] 0.99
De igual modo que para el caso anterior podemos tratar de determinar la posible asociación entre el nivel de glucosa en sangre y la HTA. Los valores medios de glucosa en sangre en hipertensos y normotensos pueden obtenerse como:
aggregate(GLUCB~HTA, data=telde, mean)
## HTA GLUCB
## 1 HTA+ 108.78086
## 2 HTA- 93.43768
Gráficamente:
a<- summarize(telde$GLUCB,by=telde$HTA, showSummary = FALSE)
EStos gráficos indican que el nivel de glucosa en sangre tiene una distribución muy asimétrica en ambos grupos. Por tanto compararlos mediante sus valores medios puede ser muy poco adecuado; además la fuerte asimetría implica que la variable bajo estudio (GLUCB) no sigue una distribución normal por lo que el test de la t de Student puede no ser válido.
En este caso resulta más conveniente realizar la comparación mediante el test de Wilcoxon-Mann-Whitney; la hipótesis nula en este test es que los valores de la variable tienden a ser similares en ambos grupos; la hipótesis alternativa es que los valores en uno de los grupos tienden a ser mayores que en el otro. Concretamente, si sospechamos que el nivel de glucosa es mayor en hipertensos, el contraste a realizar es:
wilcox.test(GLUCB~HTA, data=telde,alternative="greater")
##
## Wilcoxon rank sum test with continuity correction
##
## data: GLUCB by HTA
## W = 152680, p-value < 2.2e-16
## alternative hypothesis: true location shift is greater than 0
El p-valor indica que hay evidencia suficiente para asegurar que efectivamente en el grupo HTA+ el nivel de glucosa es sangre es mayor que en el grupo HTA-.
La siguiente base de datos:
library(foreign)
millac=read.spss("MILLAC.sav", to.data.frame=TRUE)
## Warning in read.spss("MILLAC.sav", to.data.frame = TRUE): MILLAC.sav:
## Unrecognized record type 7, subtype 18 encountered in system file
head(millac)
## ID grupo sq eLDL vLDL eTG vTG eHDL vHDL eCT vCT elpa vlpa
## 1 1 Entera-Vegetal 1 101.2 89.0 44 90 56 55 166 162 15.0 16
## 2 2 Entera-Vegetal 1 147.4 121.0 73 140 46 45 208 194 19.0 14
## 3 3 Entera-Vegetal 1 118.0 94.8 90 56 46 55 182 161 8.7 7
## 4 4 Entera-Vegetal 1 105.4 91.8 78 96 68 67 189 178 0.6 NA
## 5 5 Entera-Vegetal 1 170.4 133.0 63 60 62 57 245 202 1.4 NA
## 6 6 Entera-Vegetal 1 98.0 70.8 70 61 48 48 160 131 9.6 7
## eapoa vapoa eapob vapob
## 1 150 157 81 55
## 2 132 147 109 66
## 3 151 156 102 54
## 4 173 181 85 53
## 5 161 130 124 79
## 6 136 119 77 49
corresponde a un estudio sobre los niveles de HDL y LDL observados en niños dependiendo del tipo de leche que consumían. En concreto, el estudio se llevó a cabo mediante un diseño cruzado: los niños se dividieron aleatoriamente en dos grupos; en el primero (sq=1
) los niños consumieron leche entera durante 6 meses y leche desnatada (con grasa vegetal) durante los 6 meses siguientes; en el segundo grupo (sq=2
), durante los primeros 6 meses se consumió leche desnatada, y los siguientes 6 meses leche entera. Para cada niño se evaluaron los niveles de HDL y LDL al final de cada periodo de 6 meses. En lo que sigue denotamos por eLDL el valor de LDL tras el periodo de consumo de leche entera, y por vLDL el valor de LDL tras el consumo de leche desnatada enriquecida con grasa vegetal.
Los valores medios y desviaciones típicas de los valores de eLDL
y vLDL
observados han sido:
aggregate(eLDL~grupo,data=millac,function(x) sprintf("%.2f (%.2f)",mean(x),sd(x)))
## grupo eLDL
## 1 Entera-Vegetal 111.51 (23.50)
## 2 Vegetal-Entera 99.91 (31.09)
aggregate(vLDL~grupo,data=millac,function(x) sprintf("%.2f (%.2f)",mean(x),sd(x)))
## grupo vLDL
## 1 Entera-Vegetal 94.89 (19.95)
## 2 Vegetal-Entera 94.07 (22.83)
Así, los niveles de LDL tras la fase de consumo de leche entera (eLDL
) fueron mayores en los niños que comenzaron con leche entera (LDL medio de 111.51) que en los que comenzaron con leche desnatada (LDL medio 99.91). Asimismo, los niveles de LDL tras la fase de consumo de leche desnatada (vLDL
) fueron similares en los dos grupos (LDL medio de 94.89 entre los que hicieron la secuencia entera-vegetal, y 94.07 entre los que hicieron la secuencia vegetal-entera).
Gráficamente:
boxplot(eLDL~grupo, data=millac, main="LDL al completar la fase entera",
boxwex=0.6,col=c("orange","skyblue"))
boxplot(vLDL~grupo, data=millac, main="LDL al completar la fase vegetal",
boxwex=0.6,col=c("orange","skyblue"))
Los niños que han tomado leche en la secuencia “Entera-Vegetal” tuvieron un valor medio de LDL al terminar la fase entera de 111.51 unidades; esos mismos niños, al terminar la fase vegetal mostraron un valor medio de LDL de 94.89 unidades. ¿Es significativa esta diferencia? Es decir ¿puede explicarse por el mero efecto del azar o necesariamente debe atribuirse al efecto de haber cambiado el tipo de leche?
A diferencia de los casos anteriores para hacer esta comparación debe tenerse en cuenta que los datos están emparejados (cada niño se compara consigo mismo). Si suponemos que la variable LDL es aproximadamente normal, podemos aplicar un t.test:
eLDL=subset(millac,sq==1)$eLDL # Datos del grupo 1 (secuencia entera-vegetal) al final de la fase entera
vLDL=subset(millac,sq==1)$vLDL # Datos del grupo 1 al final de la fase vegetal
t.test(eLDL,vLDL, paired=TRUE, alternative="greater")
##
## Paired t-test
##
## data: eLDL and vLDL
## t = 5.2117, df = 46, p-value = 2.153e-06
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 11.2359 Inf
## sample estimates:
## mean of the differences
## 16.57447
lo que indica que existe evidencia suficiente para asegurar que se ha producido una disminución del LDL. La magnitud de dicha dismimución puede estimarse mediante un intervalo de confianza:
t.test(eLDL,vLDL, paired=TRUE, alternative="greater")$conf.int
## [1] 11.2359 Inf
## attr(,"conf.level")
## [1] 0.95
lo que indica que, con un 95% de confianza, que se ha producido un incremento del LDL de al menos 11.236 unidades.
Por su parte, los niños que hicieron la secuencia inversa (primero leche desnadata-vegetal y luego leche entera), tenían un nivel medio de LDL al terminar la fase vegetal de 94.07, mientras que al terminar la fase entera el LDL había subido a 99.91. Nuevamente nos hacemos la pregunta: ¿es significativa esta diferencia?
Para responderla, al igual que antes, comparamos cada niño consigo mismo a través de un t-test para muestras emparejadas:
eLDL=subset(millac,sq==2)$eLDL # Datos del grupo 2 (secuencia vegetal-entera) al final de la fase entera
vLDL=subset(millac,sq==2)$vLDL # Datos del grupo 2 al final de la fase vegetal
t.test(vLDL,eLDL, paired=TRUE, alternative="less")
##
## Paired t-test
##
## data: vLDL and eLDL
## t = -1.9392, df = 52, p-value = 0.02895
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
## -Inf -0.7803787
## sample estimates:
## mean of the differences
## -5.720755
lo que indica que efectivamente se ha producido un incremento significativo del LDL
t.test(eLDL,vLDL, paired=TRUE, alternative="less")$conf.int
## [1] -Inf 10.66113
## attr(,"conf.level")
## [1] 0.95
Por tanto podemos decir, con un 95% de confianza, que el incremento en LDL ha sido inferior a 10.66 unidades.