Veni, ViDi, Vici

Bueno, luego de realizar una pequeña introducción a las distribuciones de probabilidad, inferencia frecuentista e inferencia Bayesiana, necesarias para comprender lo que sigue, continuaré con el análisis de la salida del paquete PhViD para detección de señales en farmacovigilancia.

Teníamos nuestro espacio de trabajo cargado con el paquete PhViD y los datos de la tablita descargada con datos inventados a los fines de poder trabajar. Si no está configurado así, repetir los pasos según sugerí en el post anterior.

Ahora bien, había explicado el significado de las columnas “count”, “expected count”, “drug margin”, “event margin” y “n11/E”. Ahora nos adentraremos en la utilidad de las demás.

En primer lugar, hablaré un poco del Information Component o “IC”. Como vimos antes, la columna “n11/E” nos entrega la razón entre el valor observado y el esperado. El IC se expresa como el logaritmo natural del siguiente cociente: la probabilidad de encontrar el evento y la droga asociados, dividida por la multiplicación de la probabilidad de hallar la droga sola por la probabilidad de hallar el evento solo. Habitualmente se encuentra notada de la siguiente manera:

 Pero puede encontrarse escrito de otras maneras que tienen el mismo significado. El mismo puede calcularse también obteniendo el logaritmo natural del valor de n11/E, que es equivalente a la expresión entre paréntesis en la ecuación anterior.

En este cociente, si la ocurrencia observada en forma conjunta es igual a la esperada, el resultado será igual a 1 y su logaritmo igual a 0. Con este resultado encontraríamos que la asociación entre Droga – Evento no ocurre más frecuentemente que lo esperable y por lo tanto no se generaría una señal. Por el contrario, si es mayor a 1, tendríamos una posible señal.

El algoritmo BCPNN va más allá. Utilizando métodos de inferencia Bayesiana estima un nuevo valor para el IC, tomando como evidencia para reformularlo, los valores encontrados en la tabla. Este valor es el valor esperado del IC luego de tomar evidencia que permite estimar su valor “a posteriori”. También estima el sd “a posteriori”. Con estos valores asume (en la función por defecto, pero puede indicarse que lo haga de otra manera) que la distribución de probabilidad del IC es normal y calcula la probabilidad de obtener un valor igual o menor que 0, que sería el valor obtenido para el IC si no se genera señal alguna. Si esta probabilidad es baja se acepta la hipótesis de que hay asociación entre Droga y Evento. Esta probabilidad esta expresada en la columna postH0.

Si tenemos ganas de ver el corazón del algoritmo, podemos escribir en R:

> BCPNN

Luego de presionar “Enter” aparecerá una ristra de instrucciones que son las que se ejecutan cuando llamamos a la función. La verdad es que la mayoría es bastante críptica para alguien que no se halla acostumbrado a ver este tipo de salidas, pero hay un par de cosas que pueden notarse.

Puede verse que la función calcula dentro de su operatoria un objeto “EICb” que es la media esperada para el IC y supongo que su acrónimo significará algo como “la esperanza del IC luego de aplicar Bayes”. También, a continuación de la anterior se observa un objeto denominado “VICb” que es la varianza del IC. La raíz cuadrada de la varianza es el sd. Aquí se usa sd como equivalente del se porque estamos estimando una distribución teórica que sería equivalente a la “poblacional”, y en ese caso el sd es usado como el estimador de σ.

También, en ese listado se usa el símbolo “<-” en lugar de “=” para asignar el contenido a un objeto. Hay que saber que el primero es el nativo de R. Yo utilizo “=” por comodidad y para mayor entendimiento, pero puede pasar infrecuentemente que alguna función no arroje nada o nos de error si no utilizamos la otra forma.

Las otras columnas importantes son las denominadas “FDR” y “FNR” que son acrónimos de False Discovery Rate y False Negative Rate. Estos son conceptos introducidos por los autores del paquete y pueden verse desarrollados en los papers publicados por ellos^[1].

Básicamente, los métodos para detección de señales son muy sensibles y tienden a producir un gran número de falsos positivos, entonces se genera una variable que mide la probabilidad de obtener un falso positivo (FDR) y la de obtener un falso negativo (FNR).

Si no indicamos nada en contrario el ránking es ordenado según postH0, de menor a mayor. Esto puede cambiarse mediante el argumento “RANKSTAT”. Por defecto toma el valor igual a 1 y ordena por postH0. Acepta también el valor 2 y en este caso, produce un doble efecto: arroja una columna adicional llamada “Q_0.025(log(IC))” que indica el extremo inferior de un intervalo de confianza del 95% para el IC a posteriori, lo segundo, es que ordena el ránking por esa columna.

No obstante el orden, la regla para decidir cuáles son señales y cuáles no, es por defecto el FDR, con un punto de corte de 0.05.

Haremos ahora unas pruebas con el set de datos que viene en el paquete PhViD.

El set de datos es “PhViDdata.frame”. Lo asignaremos a un objeto que tenga un nombre un poco más manejable, por ejemplo:

> lista=PhViDdata.frame

En segundo lugar, convertiremos el objeto “lista” a un objeto en el formato que maneja el paquete:

> data=as.PhViD(lista)

Ahora usamos “BCPNN” para buscar señales, asignando la salida al objeto “Bdata”:

> Bdata=BCPNN(data)

Por último, asignaremos las señales a el objeto “signals”:

> signals=Bdata$SIGNALS

Ahora, dado que la salida será enorme el objeto signals es muy difícil de visualizar en R, y mostraré algunas funciones para poder enterarnos de lo que contiene. En primer lugar la instrucción “head”:

> head(signals)

Nos muestra el contenido de las 6 primeras filas de la tabla.

La instrucción “summary”:

> summary(signals)

Nos muestra resúmenes de las columnas, cuando son numéricos nos da estadísticos descriptivos y cuando son categóricos nos muestra las frecuencias de ocurrencias de los primeros valores.

La instrucción “str”:

> str(signals)

Nos muestra también los primeros valores pero puestos como líneas. Adicionalmente nos dice si una variable es numérica marcándola como “num” o si es categórica, marcándola como “factor”. Otra información que nos da es la dimensión de la tabla, que en este caso es de 19130 líneas x 12 columnas. Esto es bastante útil, porque el número de líneas de “signals” será igual al número de señales producidas. Si hiciéramos lo mismo con el objeto “lista” que contiene la lista de pares Droga – Evento:

> str(lista)

Obtenemos que el número total es de 102483. Vale decir, que obtuvimos:

> 19130/102483*100

[1] 18.66651

Un 18% de pares con señal positiva. Esto es lo que expresaba antes cuando decía que estos métodos daban gran cantidad de señales y que la forma de manejarlos es tomar las más intensas. Por ello, es importante el ránking ordenado de más intensa a menos intensa.

Podemos abrir una ventana con la tabla para navegarla con:

> View(signals)

Posiblemente, la forma más práctica sea exportar la tabla

> write.table(signals,"C:/PhVid/listaExp.csv",quote=F,row.names=F,sep=";")

A continuación, muestro como se genera una lista de señales dándole a “RANKSTAT” el valor de 2, para mostrar cómo se llama la función con este argumento:

> Bdata2=BCPNN(data,RANKSTAT=2)

> signals2=Bdata2$SIGNALS

Si hacemos un “str” de “signals2”, se verá que el número de señales es distinto al anterior. Esto ocurre porque el orden de la lista interviene en la creación de la variable “FDR”, que es la que determina el punto de corte.

Hay otro argumento que es importante y es “MC”. Por defecto su valor es “FALSE”. Si seteamos su valor a “MC=TRUE”, la distribución de probabilidad del IC es calculada por simulación por un método llamado Monte Carlo. Adicionalmente, el argumento “NB.MC” indica el número de simulaciones a realizar. Por defecto es 10000.

Si nuestra base es pequeña (y la mayoría los son) conviene usar este argumento para generar los parámetros calculados con menor error. Si bien el método es más lento, prioriza la precisión a la velocidad.

Finalmente, quiero decirles que no soy estadístico ni matemático ni informático. Estoy aprendiendo sobre esto casi al mismo tiempo que ustedes.

Les digo esto porque me gustaría que vuelquen sus opiniones o inquietudes en los comentarios al pie de cada post; todos serán bienvenidos. Contestaré todas las preguntas cuya respuesta sepa, y las que no, buscaremos la forma de aprenderlas.

Muchas gracias y disculpen la demora en postear. Nos vemos.

1.    Ahmed I, Thiessard F, Miremont-Salamé G, Bégaud B, Tubert-Bitter P. Pharmacovigilance Data Mining With Methods Based on False Discovery Rates: A Comparative Simulation Study. Clin Pharmacol Ther 2010 Sep;88(4):492-498.

Si es Bayes, es bueno...

Los métodos de inferencia vistos hasta el momento forman parte de lo que se da en llamar “inferencia frecuentista”. El modelo subyacente es que, verificando ciertos supuestos, es posible ajustar una distribución (normal u otra) a mis datos para inferir conclusiones acerca de la población.

La inferencia bayesiana tiene otra visión. En general parte de que la distribución “real” de la población es algo imposible de conocer. En lugar de tomar una distribución como la verdadera, se toma lo que conoce “a priori” el investigador para elegir una distribución que se ajuste a ello y trabajar inicialmente bajo ese supuesto. Luego se toma una muestra de la población y con los datos obtenidos de ella se realiza un “ajuste” de la distribución empírica, llevando a una distribución que podría ser diferente a la anterior. Es lo que se conoce como distribución “a posteriori”.

O sea que, proponemos un modelo de acuerdo a lo que conocemos, tomamos evidencia para contrastarla contra el modelo y corregimos el modelo de acuerdo a lo aprendido.

Para tener un ejemplo intentaré estimar cuál es la proporción de estudiantes de una universidad que duerme al menos 8 horas.

Este ejemplo está tomado por completo del libro “Bayesian computation with R” de Jim Albert.

Supongamos que un investigador está interesado en conocer esa proporción. El verdadero valor de la proporción le es desconocido.

De acuerdo a la información de otras publicaciones y experiencia previa, piensa 2 cosas:

·         La mejor predicción que puede hacer es de 30%

·         Es muy probable que se halle por debajo del 50%

Construimos un rango de valores posibles para la proporción. Llamaremos “p” a la proporción. En R escribimos:

> p=seq(0.05,0.95,by=0.1)

> p

 [1] 0.05 0.15 0.25 0.35 0.45 0.55 0.65 0.75 0.85 0.95

Esta instrucción construye un rango entre 0.05 (5%) y 0.95 (95%) que aumenta de a 0.1 (10%). Es un rango de valores posibles de p que nos servirá para construir una distribución a priori.

Ahora a cada uno de estos valores le otorgamos un puntaje entre 0 y 10 de acuerdo a lo que pensamos (esta escala es arbitraria). Estos puntajes son llamados “pesos”.

> pesos = c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)

> pesos

[1] 1.0 5.2 8.0 7.2 4.6 2.1 0.7 0.1 0.0 0.0

Y convertimos esos pesos a un formato de probabilidad. Para hacerlo, dividimos cada peso por la suma total de los mismos. A estas probabilidades se las llama “priors” o probabilidades “a priori”.

> prior = pesos/sum(pesos)

> prior

[1] 0.034602076 0.179930796 0.276816609 0.249134948 0.159169550 0.072664360 0.024221453 0.003460208 [9] 0.000000000 0.000000000

Esto nos permite crear una distribución que indica una probabilidad para cada valor de p. Para visualizarla podemos usar un gráfico:

> plot(p, prior, type = "h")

El gráfico obtenido nos muestra nuestra distribución “a priori”. Ahora la tarea será obtener una muestra de estudiantes. Lo hacemos y obtenemos datos de 27 estudiantes. Luego de revisar las encuestas notamos que hay 11 estudiantes que duermen más de 8 horas. Esto nos da un número de “éxitos” (11) y “fracasos” (16). Colocamos esos valores dentro de un objeto:

> data=c(11,16)

> data

[1] 11 16

El uso de la instrucción “c” es para concatenar varios valores y en este caso, asignarlos a un solo objeto.

Ahora hay que cargar un paquete de R llamado “LearnBayes”. Por supuesto, tiene que estar previamente instalado.

> library("LearnBayes")

Si no estuviese instalado nos daría un error diciendo que no encuentra el paquete. Procedemos como en el post que explica sobre instalación de paquetes de R y lo instalamos. Luego volvemos a tipear la instrucción anterior.

Luego de cargado el paquete escribimos el siguiente comando:

> post = pdisc(p, prior, data)

> post

[1] 1.451714e-08 2.256022e-03 1.291349e-01 4.767910e-01 3.338350e-01 5.587820e-02 2.098297e-03 [8] 6.642742e-06 0.000000e+00 0.000000e+00

Esto nos da las probabilidades a posteriori. Es bastante claro que son diferentes a las probabilidades a priori.

Para visualizarlas en forma comparativa escribimos una serie de instrucciones que no explicaré en detalle, pero nos permiten ver una gráfica sobre la otra. En caso de que a alguien le interese la explicación detallada me lo puede pedir. No obstante, se puede obtener ayuda de R escribiendo “?” seguido de la instrucción, o help(“instrucción”).

> library(lattice)

> PRIOR=data.frame("prior",p,prior)

> POST=data.frame("posterior",p,post)

> names(PRIOR)=c("Type","P","Probability")

> names(POST)=c("Type","P","Probability")

> data=rbind(PRIOR,POST)

>xyplot(Probability~P|Type,data=data,layout=c(1,2),type="h",lwd=3,col="black")

Esto nos muestra que las 2 distribuciones, la a priori construida con las suposiciones del investigador y la a posteriori, calculada usando la a priori y corrigiendo de acuerdo a las observaciones, son diferentes.

Existen funciones matemáticas para calcular las distribuciones a posteriori a partir de las a priori. Estas funciones son diferentes en cada situación y mi objetivo en este post no es ser extenso en éste tema.

Baste decir, que con esta metodología es posible estimar distribuciones pero también parámetros, como medias, sd y otros.

Con estos conocimientos podemos ir, ahora sí, a continuar con el aprendizaje del paquete PhViD para detección de señales.

Yo infiero, tu infieres

Luego de repasar algunos rudimentos de distribuciones estadísticas, es conveniente ver algo de inferencia estadística. Es la parte de la estadística que permite deducir propiedades de una población, a partir de una muestra de la misma.

Acá hay algunos conceptos que habitualmente conducen a error. Básicamente, como decíamos antes, cuando estudiamos una muestra una de las cosas que hay que realizar es una serie de pruebas para establecer si la misma proviene de una población donde la variable tiene una distribución normal. Esto implica que existe una población teórica (e imposible de aprehender) desde la cual se toma la muestra. En esta muestra verificaremos propiedades similares a la de la población, pero no idénticas. Por ejemplo, si calculamos la media en la muestra de glucemias con la que trabajamos anteriormente veremos que:

> mean(glucemias)

[1] 90.03842

Su promedio no es exactamente 90, que es lo esperado conociendo la media real, pero está muy cerca. Esto nos puede llevar a una pregunta. Si tomamos muchas muestras de la misma población y calculamos sus promedios, ¿cuánto se alejarán éstos de la media poblacional? Es aquí que surge la idea conceptualmente diferente de la que hablaba antes. Existe una distribución que se puede ajustar a los valores de nuestra muestra para describirla, pero lo que nos interesa en realidad, es una distribución teórica que muestra cómo se distribuirían los promedios si tomáramos muchas muestras respecto de la media poblacional. Inicialmente, es fácil imaginar que la mayoría de los promedios estarán cerca de la media y si bien obtendremos por azar valores alejados, éstos se irán haciendo menos y menos probables a medida que nos alejamos de ella. Esto es muy importante, existe una distribución teórica de medias muestrales alrededor de la media poblacional. Esta distribución tiene a su vez una media (la media poblacional o “verdadera” media) y un desvío estándar, que en general se denotan usando las letras griegas µ y σ para diferenciarlas y resaltar que en general son valores teóricos desconocidos y a estimar.

Éste es el razonamiento básico subyacente en lo que solemos ver como Intervalo de Confianza (ojo que su acrónimo en español “IC” puede confundirse con el de Information Component que nos interesa, pero no es lo mismo).

Si tuviéramos información acerca de la media poblacional y su desvío estándar, podríamos calcular que tan probable es obtener un promedio determinado para una muestra tomada. Como no lo tenemos, lo inferimos. Quiero decir, vamos al revés. En lugar de ir de la información de la población para evaluar la muestra, vamos de la información que hay en la muestra para evaluar cómo será la población.

Si tuviéramos los datos de la población, calcularíamos la media y el desvío estándar de ella, mediríamos a qué distancia en términos de sd se encuentra la media de nuestra muestra y buscaríamos la probabilidad de encontrar un valor así de alejado de la media poblacional. Como no conocemos la media y el sd poblacional, usamos la media muestral para estimar la poblacional y para estimar el sd poblacional hacemos una corrección del sd muestral. Esta corrección hará que el sd sea más estrecho, como sería de esperar si tuviéramos toda la población y obviamente, está relacionada al tamaño de la muestra (n).

 Este nuevo valor del sd es llamado Error Standard (se) y se trata del valor que usamos para estimar el sd poblacional. Entonces con esto podemos calcular un intervalo de confianza para una muestra. SI luego del análisis preliminar concluimos que nuestra muestra proviene de una población con distribución de valores normal, podemos usar aproximadamente 2 errores estándar (sabemos por lo visto en el post anterior que el valor exacto es de 1.959964 sd) hacia la derecha y 2 errores estándar hacia la izquierda del promedio para establecer un intervalo de confianza de 95%. Entre las características que podemos usar para verificar si se cumple el principio de normalidad, uno muy importante es el tamaño de la muestra. Cuanto mayor el tamaño muestral, menos dependemos del resto de los supuestos.

Este intervalo de confianza tiene el siguiente significado: si tomáramos muchas muestras con el mismo método que tomamos la primera, el 95% de las medias se encontraría en ese rango y además, lo que es más importante, en ese rango se encuentra la media poblacional con un 95% de probabilidad.

Veamos esto con las glucemias. Usaremos la función “mean” y “length” ya vistas, la función “sd” que calcula el desvío estándar y la función “sqrt” que nos da la raíz cuadrada.

> nG=length(glucemias)

> meanG=mean(glucemias)

> sdG=sd(glucemias)

> seG=sdG/sqrt(nG)

> sdG

[1] 9.519673

> seG

[1] 0.9519673

Vemos que el sd es de 9.52, cercano a 10 que es el valor poblacional, pero el se es mucho menor, cercano a 0.92. Para obtener el intervalo de confianza acerca de nuestra media hacemos lo siguiente:

> LimInf=meanG-1.959964*seG

> LimSup=meanG+1.959964*seG

> LimInf

[1] 88.1726

> LimSup

[1] 91.90424

El límite inferior es de 88.17 y el superior de 91.90. Como dijimos antes, esto quiere decir entre otras cosas que la media poblacional “real” se encuentra en ese rango con un 95% de probabilidad. Sabemos que esto es lo que ocurre en la realidad, pues la media poblacional es de 90 mg/dl, según lo predeterminado por la simulación con la que generamos los datos muestrales.

También, nos dice otra cosa muy importante: Si en otro momento tomamos otra muestra y al calcular el promedio encontramos que se encuentra fuera de ese rango, podemos tener la “confianza” de que es improbable que esos valores provengan de la misma población. Mejor dicho, la probabilidad de que provengan de la misma población es inferior al 5%.

Bueno, hasta aquí los conceptos de estadística “clásica” que son importantes para entender lo que viene de generación de señales con el paquete PhViD.

No obstante, para entenderlo más acabadamente hay que tener nociones de inferencia Bayesiana.

Se me caen las medias!

Voy a hablar un poco acerca de distribuciones de probabilidad. Este es un tema extenso así que solamente trataré los puntos que me parecen necesarios para entender los demás.

Como vimos en el post sobre nociones de independencia, en las tablas de 2 x 2 existe un valor observado para cada celda interna, que es el que anotamos a partir de las observaciones, y un valor esperado que podemos calcular y es a su vez el valor más probable bajo la hipótesis de que las variables son independientes.

Entonces, queda configurada una “distancia” entre el valor esperado y el observado. Intuitivamente puede verse que a mayor diferencia del valor observado respecto del esperado, menos probable se hace el mismo si las variables fueran independientes y por lo tanto se hace cada vez más difícil aceptar esta hipótesis. Cuando la distancia es demasiada (y para ello se define un punto de corte), rechazamos la hipótesis de independencia y aceptamos que hay asociación entre las variables, pues de otro modo sería sumamente improbable obtener esa configuración de celdas.

Ahora bien: ¿Cómo se mide esa distancia y cómo se hace para saber el valor de probabilidad de obtener esa diferencia? Aquí entran las distribuciones de probabilidad. Usándolas, podemos traducir entre distancia medida en términos comunes y probabilidad asociada a la misma.

Para tablas de 2 x 2, en general se usa la distribución de chi cuadrado, porque las propiedades que tienen estas tablas establecen que este método sea adecuado la mayor parte de las veces. Pero otras veces, lo adecuado es elegir otras distribuciones para medir la distancia observada y eso depende de la situación a la que nos enfrentemos.

Para dar una noción más intuitiva del uso de las distribuciones de probabilidad, voy a ejemplificar con la distribución llamada Normal que es posiblemente la más escuchada.

Es de esperar que la mayoría de nosotros haya visto alguna vez una de esas reglas que sirven para medir en diferentes escalas. Por ejemplo, una regla que de un lado tiene una escala en centímetros y del otro en pulgadas.

Una regla así, nos podría servir como método para traducir entre una escala y otra. Es fácil ver en la imagen que si medimos algo que mide 1 pulgada, medirá aproximadamente 2,5 centímetros. Algo así pasa con las distribuciones de probabilidad, vale decir que podríamos medir una distancia en centímetros y traducir a probabilidades. Para eso, vamos a usar una herramienta que se consigue en cualquier ferretería estadística. ;)

Esta regla tendrá, al menos inicialmente, una escala en centímetros. Pero posee algunas características particulares. En primer lugar, el valor 0 no está colocado en un extremo sino en el centro de la misma. En segundo lugar, los centímetros son positivos hacia la derecha del cero pero se hacen cada vez más negativos hacia la izquierda del mismo.

Como esta no es una regla cualquiera, tiene además la particularidad de que se le pueden acoplar diferentes escalas en el borde opuesto al de los números, de forma de poder traducir de escala en centímetros a cualquier otra.

En esta oportunidad, le acoplaremos una curva de probabilidad Normal:

Como ven, esta nueva escala es bastante rara. Básicamente, a diferencia de lo que ocurre con la escala en centímetros, su grosor es máximo en el centro o valor 0 y éste se va adelgazando a medida que nos alejamos del mismo hasta hacerse mínimo más allá de los 3 centímetros de distancia respecto del centro. Obviamente, el área cubierta por la figura en cada segmento también disminuye. También podemos ver que la figura es simétrica hacia ambos lados.
Lo interesante (como veremos) es que si asumimos que el área total de la figura es 100%, su forma es tal que es posible saber exactamente qué porcentaje del área queda encerrada en cada segmento medido en centímetros. En primer lugar, es fácil darse cuenta que cualquiera de las 2 mitades contiene un área de 50%. Entre 0 y 1 cm queda determinada un área de 34% del total y dada su simetría lo mismo ocurre entre 0 y -1 cm. Por lo tanto, entre -1 y 1 cm queda un área de 68%.

Varios rangos que podríamos verificar serían:

Rango en cm	Porcentaje del área
0 a 1	34%
-1 a 1	68%
0 a 2	47,7%
-2 a 2	95,4%
0 a 3	49,9%
-3 a 3	99,7%
-1 a 3	84%
Todo el rango menor o igual a cero	50 %
Todo el rango menor o igual a 2	97,7%

La verdad es que si bien estuvimos hablando en términos de porcentajes, también estuvimos hablando en términos de probabilidades. El principio fundamental es que no debemos nunca usar un destornillador para sacar un clavo. Vale decir, cada herramienta tiene una situación de aplicación para lo cual fue diseñada.

Así, lo más importante antes de hacer una medición con la herramienta que mostramos, es saber si nuestros datos pueden ser medidos con la misma o debe usarse otra. Sabemos por lo que leímos en muchas publicaciones, que muchos valores en medicina tendrían una gráfica de distribución normal y por lo tanto simétrica alrededor de la media. (La palabra “tendría” es condicional aquí, pues hay situaciones en las que esto no se cumple).

Vamos a usar R para simular datos con distribución normal y reforzar lo que venimos hablando.  Pretendamos que estamos midiendo valores de glucemia en una población sana en mg/dl y que el rango normal va de 70 a 110 mg/dl. Usaremos el valor intermedio del rango como media (90 mg/dl) y un desvío estándar de 10 mg/dl.

La función “rnorm” sirve para generar números aleatorios según la distribución normal y toma como argumentos la media y el desvío estándar (sd), además del número de observaciones o “casos” que queremos. Pediremos 100 casos:

> n=100

> media=90

> sd=10

> glucemias=rnorm(n,media,sd)

Luego, graficaremos la densidad de distribución de nuestros datos de la siguiente manera:

plot(density(glucemias))

Obtendremos una gráfica más o menos así:

Posiblemente la gráfica sea medio torcida, pero vemos que es bastante parecida a la normal (si en lugar de un n de 100 usáramos 1000 o 10000, veríamos que cada vez es más parecida).

Como suponemos que es aceptable usar la distribución normal para estudiar los datos (aquí este supuesto es obviamente correcto pues los datos fueron generados por simulación) mediremos la probabilidad de obtener valores dentro de determinados rangos y los compararemos con la tabla de probabilidad normal aportada más arriba.

Para obtener los valores dentro de un determinado rango, tomaremos nuestro objeto “glucemias” que contiene nuestros datos y lo “cortaremos” según necesidad. Para tomar segmentos de una lista de valores contenida en un objeto, se usa el nombre del objeto seguido de corchetes (“[ ]”), y dentro de los corchetes irá la instrucción que indica dónde cortaremos.

> glucemias[glucemias<=90]

Esto nos dará la lista de valores de glucemia que tienen valor menor o igual a nuestra media, que es el valor de 90 mg/dl. Como tenemos 100 casos totales, el número de casos en un segmento será también el porcentaje del total de casos contenido en ese segmento. Para saber el número de casos que hay en la lista de valores sin tener que contarlos a mano usamos la instrucción “length”, envolviendo a la instrucción anterior de la siguiente manera:

> length(glucemias[glucemias<=90])

Rango en sd	Porcentaje del área esperada	Cantidad obtenida en el rango	Instrucción
0 a 1	34%	31	length(glucemias[glucemias<=100])- length(glucemias[glucemias<=90])
-1 a 1	68%	64	length(glucemias[glucemias<=100])- length(glucemias[glucemias<=80])
0 a 2	47,7%	49	length(glucemias[glucemias<=110])- length(glucemias[glucemias<=90])
-2 a 2	95,4%	99	length(glucemias[glucemias<=110])- length(glucemias[glucemias<=70])
0 a 3	49,9%	49	length(glucemias[glucemias<=120])- length(glucemias[glucemias<=90])
-3 a 3	99,7%	100	length(glucemias[glucemias<=120])- length(glucemias[glucemias<=60])
-1 a 3	84%	82	length(glucemias[glucemias<=120])- length(glucemias[glucemias<=80])
Todo el rango menor o igual a la media	50 %	51	length(glucemias[glucemias<=90])
Todo el rango menor o igual a 2 sd	97,7%	100	length(glucemias[glucemias<=110])

Es interesante verificar que al obtener una muestra de una población, a pesar de tener la variable medida una distribución expresamente normal en esa población, su distribución en la muestra será aproximada pero no exactamente normal. Habitualmente el procedimiento será inverso, se tomará una muestra de una población y deberán realizarse algunas evaluaciones para determinar si se cumplen los supuestos que indican que la distribución normal es pasible de ser aplicada sin incurrir en violaciones de los mismos que, en general, conducirían a resultados y conclusiones erróneas.

¿Pero, para qué utilicé la analogía de los centímetros? Y ¿Por qué en una tabla puse los rangos en centímetros y en otra en sd (desvíos estándar)? Existe una distribución normal especial, denominada “Normal(0,1)”, en donde la media es igual a cero y el desvío estándar igual a 1. Si bien cualquier probabilidad puede ser calculada computacionalmente sin recurrir al uso de tablas, es útil saber que las tablas de probabilidad normal que a veces se ven publicadas y algunos tuvimos que usar alguna vez, están expresadas según la Normal(0,1). Aunque tengamos una media y sd diferentes, cualquier valor puede ser fácilmente convertido a esa escala mediante una cuenta en donde restamos el valor de la media y luego lo dividimos por el desvío estándar:

Así, sabemos a cuántos desvíos estándar está nuestro valor, respecto de la media y podemos buscarlo en las tablas. En R podemos usar el comando “pnorm” con el resultado de ese valor, para obtener la probabilidad de obtenerlo. El resultado será evaluado como si se tratara de uno obtenido de una Normal(0,1) (salvo que indiquemos otros argumentos). Por ejemplo, si la cuenta anterior diera un resultado de 2 (sería el caso si con las glucemias usáramos el valor de 110 mg/dl), sabríamos que este valor está a 2 sd respecto de la media. Usamos “pnorm” para establecer la probabilidad de obtener un valor así, tomado de una distribución normal:

> pnorm(2)

[1] 0.9772499

Tal vez hubiéramos esperado que el valor estuviera cerca del 95% (a esta altura la mayoría se imaginará que cuando hablo de 95% y 0.95 me refiero a lo mismo. En realidad las probabilidades se expresan en valores de 0 a 1. Para transformarlas en porcentaje se multiplican por 100), de acuerdo a los conocimientos que suelen transmitirnos. En realidad, la función da la probabilidad de obtener un valor igual o menor a 2, no entre ±2sd o lo que es lo mismo, entre 2 y -2. Si queremos obtener esto último, debemos restarle la probabilidad de obtener un valor menor o igual a -2, de la siguiente manera:

> pnorm(2)-pnorm(-2)

[1] 0.9544997

Como ven, está cercana a 95%, aunque no exacta. Si quisiéramos obtener el valor para ±sd que nos da exactamente el 95% de probabilidad, usamos “qnorm” indicándole el valor 0.975 (o sea 97,5%, que es el valor que deja 2,5% hacia la derecha, pues queremos que nos quede un 2,5% a la izquierda).

> qnorm(0.975)

[1] 1.959964

Ahora volvemos a utilizar “pnorm” con ese valor:

> pnorm(1.959964)-pnorm(-1.959964)

[1] 0.95

Exactamente 95%.

Bueno, acá quiero que respiren un poco, se estiren y despejen un poco la cabeza. Vamos a pasar a un punto importante.

Cuando hacemos cualquier estudio, la mayoría de las veces queremos generalizar el resultado. Vale decir, puede ser que las conclusiones nos interesen solamente para ser aplicadas para esos pacientes y ese momento particular, pero habitualmente queremos que sean útiles para usarlas en otros pacientes. Esto implica que las conclusiones serán utilizables para aplicar a otros pacientes similares en otras regiones, en bases de datos históricas, y en pacientes que aún no han llegado, pero llegarán en el futuro. Esto quiere decir que cuando tomo una muestra, no solamente lo hago porque es impracticable estudiar a toda la población de interés por una cuestión de recursos sino que esa población es una población teórica que la mayoría de las veces no existe al momento de realizar el estudio.

Hasta el momento, con lo que tenemos, solamente podemos describir nuestra muestra en términos de media y desvío estándar. Ahora, si queremos evaluar cuán aplicable es mi resultado a otros pacientes no incluidos en el estudio, debemos usar otro razonamiento y realizar un salto conceptual para entrar en el terreno de la inferencia estadística.