6.6 Aplicaciones

La creación de la tabla de datos se puede hacer como se indica a continuación, a partir de una tabla creada con la matriz de comandos() :

> tab < – as.tabla (matriz (c (815,115,208,327), nrow = 2, byrow = TRUE,

dimnames = lista (EST = c(«+»,»−»),

CAD = c(«+»,»-«))))

> tab

Cabe señalar que la tabla se reorganizó para corresponder a las notaciones habituales, con eventos «positivos» en la primera fila (típicamente, la exposición) y la primera columna (típicamente, la enfermedad) de la tabla.

Entonces sería posible calcular todas las cantidades requeridas usando muy pocas operaciones aritméticas con R. Sin embargo, el paquete epiR contiene todas las herramientas necesarias para responder a las preguntas epidemiológicas en las tablas 2 × 2. Por lo tanto, el comando epi.tests() proporciona valores de prevalencia,sensibilidad / especificidad, así como valores predictivos positivos y negativos:

> biblioteca(epiR)

> epi.pruebas (tab)

Con respecto al área bajo la curva, también se puede emplear un comando externo, roc.de.tabla(), en el paquete epicalc:

> library(epicalc)

> roc.de.tabla (tab, gráfico = FALSO)

Importar los datos de idh.dat se logra a continuación:

> bp < leer.tabla(«hdis.dat», header = TRUE)

> str (bp)

Como se puede comprobar, no hay etiqueta asociada a los niveles de las variables de interés (bpress para la presión arterial, chol para la tasa de colesterol). Para generar y asociar etiquetas, se pueden insertar los siguientes comandos:

> blab < – c(«<117”,»117-126”,»127-136”,»137-146”,

«147-156”,»157-166”,»167-186”,»>186”)

> clab < – c(«<200”,»200-209”,»210-219”,»220-244”,

«245-259”,»260-284”,»>284”)

> bp$bpress < – factor(bp$bpress, labels = blab)

> bp$chol < – factor (bp cho chol, labels = clab)

El último comando convierte las variables originales en variables cualitativas y asocia a sus modalidades las etiquetas definidas por blab y clab. Para verificar que la base de datos está ahora en la forma deseada, los siguientes comandos son útiles:

> str (bp)

> summary(bp)

La tabla de las frecuencias relativas expuestas en la instrucción ahora se puede reproducir:

> round (xtabs (hdis / total bpress + chol, data = bp), 2)

Dado que nos vamos a centrar en la relación entre la presión arterial y el infarto de miocardio solamente, es necesario agregar datos sobre los niveles de colesterol. En otras palabras, es necesario acumular los valores para cada nivel de presión arterial, independientemente de los niveles de colesterol. También se debe aprovechar la oportunidad para cambiar el nombre de los niveles de la variable bpress utilizando los centros de los intervalos de clase:

> blab2 < – c(111.5,121.5,131.5,141.5,151.5,161.5,176.5,191.5)

> bp$bpress < rep(blab2, cada uno = 7)

> dfrm < – agregado(bp,

lista(bpress = bp), suma)

es el último comando aggregate(), que hace que sea posible combinar los datos: suma todos los valores (sum) de la variable chol que no están incluidos en la lista de variables que serán objeto de análisis. Al mismo tiempo, los resultados se almacenan en una nueva base de datos llamada dfrm. Tenga en cuenta que en este caso no hacemos uso de la notación de fórmulas, sino de la sintaxis alternativa para aggregate() (variable de respuesta seguida de una lista que consta de una o más variables explicativas). Se puede obtener una visión general de los datos agregados dando head() como entrada:

> head(dfrm, 5)

Cuando una proporción p está disponible, su valor en una escala cuyas unidades son logits viene dado por el registro de relación(p/(1 − p)). De ahí se pueden derivar los siguientes comandos para convertir las proporciones, calculadas como hdis / total en unidades logit:

> logit < – function(x) log(x/(1-x))

> dfrm$prop < dfrm$idh/dfrm$total

> dfrm$logit < – logit(dfrm$idh/dfrm$total)

se observa que una pequeña función ha sido definida por la conversión de los valores de x, que aquí se supone que es un proporciones, en su equivalente log(x/(1 − x)). Igualmente, uno podría escribir:

> log((dfrm$idh/dfrm$total)/(1-dfrm$idh/dfrm$total))

Los valores calculados se han añadido directamente al marco de datos bajo el nombre de pilar y logit.

El modelo de regresión logística está escrito de la siguiente manera:

> glm(cbind(hdis,total-hdis) bpress,data = dfrm,family = binomial)

La formulación utilizada, cbind(hdis, total-hdis) bpress, tiene en cuenta el hecho de que estamos accediendo a datos agrupados y no a respuestas individuales. El comando glm () con la opción family = binomial corresponde a una regresión logística, que, sin entrar en detalles, utiliza por defecto la función logit como enlace canónico.

Los coeficientes de regresión se pueden mostrar utilizando summary ():

> summary (glm (cbind (hdis, total-hdis) bpress, data = dfrm,

family = binomial))

El resultado anterior incluye la información esencial para responder a la pregunta de la significación estadística de la asociación entre la presión arterial y la probabilidad de ataque cardíaco: el coeficiente de regresión (en la escala de probabilidades logarítmicas) es igual a 0,024 y es significativo en el umbral habitual del 5% (véase la columna Pr(>|z|)).

La probabilidad de sufrir un ataque cardíaco de acuerdo con los diferentes niveles de presión arterial considerados se obtiene de la siguiente manera:

> glm.res < – glm (cbind(hdis, total-hdis) bpress, data = dfrm,

family = binomial)

> predict (glm.res)

Cabe señalar que los resultados intermedios generados por R se han almacenado con el nombre glm.res antes de usar el comando predict (). Las predicciones generadas por R se expresan en forma de logit y es posible comparar los logits observados y predichos:

> cbind(dfrm, logit.predit = predecir (glm.(res))

Prácticamente todos los elementos están disponibles para representar gráficamente las proporciones de infarto de miocardio observado y predicho de acuerdo con el nivel de presión arterial. Faltan las predicciones del modelo expresadas en forma de proporciones y no en la escala logodds. Por lo tanto, puede ser deseable dibujar la curva de predicción, es decir, la probabilidad de tener un ataque cardíaco de acuerdo con la presión arterial, sin limitar esta última a los ocho valores observados para la variable bpress. Aquí sigue una posible solución:

> dfrm$prop.predit < – predict (glm.res, type = «respuesta»)

> f < – function(x)

1/(1 + exp(-(coef(glm.res) + coef(glm.res)*x)))

> xyplot(hdis/total ˜ bpress, data = dfrm, aspect = 1.2, cex = .8,

xlab = «Blood pressure”,ylab = «Infarction likelihood”,

panel = function(x, y, …) {

panel.xyplot(x, y, col = «gray30”, pch = 19, …)

panel.curve(f, lty = 3, col = «gray70”)

panel.points(x, dfrm$prop.predit, col = «gray70»,…)

})

Dado que los datos brutos (individuales) no están disponibles, es necesario trabajar directamente con la tabla de frecuencias proporcionada en la declaración del problema:

> alcohol < – matrix(c(666,104,109,96), nr = 2,

dimnames = list(c(«Control”,»Case”),

c(«<> = 80”)))

> alcohol

Regarding the odds ratio, the command oddsratio() from the package vcd will be used:

> biblioteca(vcd)

> oddsratio (alcohol, log = FALSE)

La opción log = FALSE garantiza que el resultado devuelto corresponda correctamente a una relación de probabilidades y no al registro de la relación de probabilidades. Para obtener un intervalo de confianza asintótico, se empleará lo siguiente:

> confint(oddsratio(alcohol, log = FALSE))

De manera similar, el resumen(oddsratio (alcohol)) podría usarse para realizar una prueba de hipótesis sobre la relación de probabilidades logarítmicas(H0 : log (θ) = 0).

La prop de comando.test () se puede utilizar para probar la hipótesis de que la proporción de personas con un consumo diario es ≥ 80 g y es idéntica tanto en los casos como en los controles, indicando los valores observados a partir de la tabulación cruzada dada en la declaración del problema:

> prop.prueba(c(96,109), c (200,775), correct = FALSE)

Esta prueba es exactamente equivalente a una prueba Z para probar la diferencia entre dos proporciones estimadas a partir de los datos (si no se está empleando una corrección de continuidad).

Con respecto al modelo de regresión logística, la tabla de contingencia debe transformarse en una tabla de datos en la que las dos variables cualitativas (enfermedad y exposición) estén representadas adecuadamente, es decir, con un marco de datos. Esto se puede hacer usando el comando melt() del paquete reshape2 que hace posible transformar los datos en forma tabular en formato «largo». En esta etapa, los niveles de la variable enfermedad deben recodificarse con 0 (control) y 1 (casos), aunque esto no es realmente necesario y el nivel 0 debe considerarse como la categoría de referencia (que facilita la interpretación de los resultados):

> biblioteca(reshape2)

> alcohol.df < – melt(alcohol)

> nombres(el alcohol.df) < – c(«enfermedad», «exposición», «n»)

> niveles(el alcohol.df disease enfermedad) < – c(1,0)

> alcohol.df disease disease < – relevel (alcohol.df disease enfermedad, «0»)

El modelo de regresión logística se escribe de la siguiente manera:

> glm.res < – glm (exposición a enfermedades, datos = alcohol.df,

family = binomial, pesos = n)

> resumen(glm.res)

El resultado de interés es la fila asociada a la exposición > = 80, ya que proporciona información sobre el valor del coeficiente de regresión asociado a la variable de exposición estimada por R, con su error estándar, así como el valor del estadístico de prueba. En este caso, el coeficiente de regresión se interpreta como el logaritmo de la razón de probabilidades. Tenga en cuenta que se obtendrían exactamente los mismos resultados intercambiando el papel de las variables en la formulación anterior, enfermedad por exposición.

La razón de probabilidades calculada anteriormente se puede encontrar a partir del coeficiente de regresión asociado al factor de interés (exposición), así como su intervalo de confianza del 95%, de la siguiente manera:

> exp(coef(glm.res))

> exp (confint (glm.res))

La segunda solución consiste en considerar el número de casos y el número total de individuos:

> alcohol2 < datos.marco (expos = c («<> = 80»), case = c(104,96),

total = c(104 + 666, 96 + 109))

> resumen (exposiciones glm(cbind(caso, caso total), datos = alcohol2,

familia = binomio))

4) Estos datos provienen de un estudio que busca establecer la validez pronóstica de la concentración de creatincinasa en el cuerpo sobre la prevención de la aparición de infarto de miocardio .

Dado que el nombre de las variables no está en el archivo de datos, será necesario asignarlas inmediatamente después de importar los datos:

> sck < – lectura.tabla(«sck.dat», header = FALSE)

> nombres(sck) < – c(«ck», «presión», «abs»)

> resumen(sck)

El número total de sujetos que corresponde a la suma de las frecuencias para las dos variables pres y abs, que es:

> sum(sck)

o, equivalentemente,: suma(sck$pres) + sum(sck$abs) (pero no sck$pres + sck$abs).

Para calcular las frecuencias relativas de estas dos variables, es necesario conocer las cantidades totales por variable. Estos se pueden obtener utilizando el comando apply y operando por columnas:

>ni < – apply(sck, 2, sum)

En base a esto, es suficiente dividir cada valor de las variables pres y abs por los valores calculados previamente. Los valores obtenidos se almacenarán en dos nuevas variables en la misma tabla de datos:

> sck$pres.la proposición < – sck$pres/ni

> sck$abs.prop < – sck abs abs/ni

Es posible verificar que los cálculos son correctos: la suma de los valores para cada variable ahora debe ser igual a 1:

> apply(sck, 2, sum)

Entonces las proporciones obtenidas se representan simplemente en un solo gráfico, teniendo en cuenta los valores variable ck como abscisas:

> xyplot(pres.prop + abs.prop ck, data = sck, type = c («b»,»g»),

auto.key = TRUE, ylab = «Frequency»)

La instrucción type = c («b»,» g») significa que queremos que los puntos se muestren conectados por líneas («b» = » o » + «l’), así como una cuadrícula («g»).

El modelo de regresión para los datos agrupados se escribe como:

> glm.res < – glm(cbind(pres,abs) ck, data = sck, family = binomial)

> resumen(glm.res)

Las predicciones, expresadas en forma de probabilidades y no en la escala log-odds, se obtienen mediante el comando predecir especificando la opción type = «response»de la siguiente manera:

> glm.pred < – predict (glm.res, type = «response»)

> nombres (glm.pred) < – sck c ck

Considerando que las probabilidades ≥ 0,5 designan individuos enfermos, se obtiene la siguiente distribución:

> glm.pred

Se concluye que las personas serán consideradas como enfermas, según este modelo, para valores ck de 80 o más. Esto se verifica fácilmente con un gráfico donde se representan las probabilidades predichas basadas en los valores de la variable ck:

> sck sick sick < – sck pre pres/(sck pre pres + sck abs abs)

> xyplot(glm.pred ~ sck$ck, type = «l»,

ylab = «Probabilidad», xlab = «ck»,

panel = function(…) {

panel.panel xyplot (panel)

.xyplot (sck c ck, sck sick sick, pch = 19, col = «grey»)

})

En un primer paso, es necesario «descomprimir» los datos agrupados y crear una tabla con dos columnas: la primera representando la variable ck y la segunda representando la presencia o ausencia de enfermedad. Usaremos los recuentos por subgrupo calculados previamente y disponibles en la variable ni:

> sck.expanda < datos.marco(ck = c(rep(sck$ck, sck$pres),

rep(sck$ck, sck$abs)),

enfermo = c(rep(1,ni), rep(0,ni)))

> tabla(sck.expand sick sick)

> con(sck.expand, tapply(sick, ck, sum))

Los dos últimos comandos están destinados a garantizar que se encuentren los mismos recuentos y que la distribución de personas enfermas de acuerdo con el valor de ck sea correcta. También se puede verificar que se obtienen los mismos resultados con respecto a la regresión logística y al mismo tiempo se pueden agregar los valores predichos en la tabla de datos anterior (se podría seguir el mismo procedimiento que anteriormente y se podrían replicar las predicciones, pero esto es más simple de esta manera):

> glm.res2 < – glm (sick ck, data = sck.expand, family = binomial)

> sck.expand prediction prediction < – ifelse (predict (glm.res2,

type = «response») > = 0.5,

1, 0)

> con(sck.expand, table (sick, prediction))

El último comando permite mostrar una matriz de confusión en la que se cruzan los diagnósticos reales y los predichos por el modelo de regresión. Las tasas de clasificación correctas se pueden comparar cuando variamos el umbral de referencia:

> classif.tab < – con (sck.expand, table(sick, prediction))

> sum(diag (classif.tab)) / sum (classif.para mostrar la curva ROC, haremos uso del paquete ROCR:

> library(ROCR)

> pred < – prediction (predict (glm.res2, type = «response»),

sck.expand$sick)

> perf < – performance(pred, «tpr”,»fpr”)

> plot(perf, ylab = «Sensitivity”, xlab = «1-Specificity”)

> grid()

> abline(0, 1, lty = 2)

The value of the area under the curve is obtained as follows:

> performance(pred, «auc”)@»y.values”