Estadística archivos - Página 2 de 9 - Ciencia sin seso... locura doble

Tres patas de un gato

Publicado el 28 febrero, 2017 | por Manuel Molina | Deja un Comentario

Lo de buscarle tres pies al gato, o tres patas, es un dicho muy popular. Parece que se dice que busca tres pies a un gato aquél que trata de demostrar alguna cosa imposible, generalmente con tretas y engaños. En realidad, el refrán inicial hacía referencia a buscar cinco pies en lugar de tres. Esto parece más lógico, ya que como los gatos tienen cuatro patas, encontrarles tres de ellas es cosa fácil, pero encontrar cinco es algo imposible, a no ser que consideremos la cola del gato como otro pie, lo cual no tiene mucho sentido.

Pero hoy no vamos a hablar de gatos con tres, cuatro o cinco pies. Vamos a hablar sobre algo un poco más etéreo, como son los modelos multivariables de regresión lineal múltiple. Este sí que es un gato con multitud de pies, pero nosotros nos vamos a fijar únicamente en tres de ellos que reciben los nombres de colinealidad, tolerancia y factor de inflación (o incremento) de la varianza. Que nadie se desanime, es más fácil de lo que puede parecer de entrada.

Ya vimos en una entrada anterior cómo los modelos de regresión lineal simple relacionaban dos variables entre sí, de forma que las variaciones de una de ellas (la variable independiente o predictora) podían servir para calcular cómo iba a variar la otra variable (la variable dependiente). Estos modelos se representaban según la ecuación y = a + bx, donde x es la variable independiente e y la dependiente.

Pues bien, la regresión lineal múltiple añade más variables independientes, de tal manera que permite hacer predicciones de la variable dependiente según los valores de las variables predictoras o independientes. La fórmula genérica sería la siguiente:

y = a + bx₁ + cx₂ + dx₃ + … + nx_n, siendo n el número de variables independientes.

Una de las condiciones para que el modelo de regresión lineal múltiple funcione adecuadamente es que las variables independientes sean realmente independientes y no estén correlacionadas entre sí.

Imaginad un ejemplo absurdo en el que metemos en el modelo el peso en kilogramos y el peso en libras. Ambas variables variarán del mismo modo. De hecho el coeficiente de correlación, R, será 1, ya que prácticamente las dos representan la misma variable. Ejemplos tan tontos es difícil verlos en los trabajos científicos, pero hay otros menos evidentes (como incluir, por ejemplo la talla y el índice de masa corporal, que se calcula a partir del peso y de la talla) y otros que no son evidentes en absoluto para el investigador. Esto es lo que se llama colinealidad, que no es más que la existencia de una asociación lineal entre el conjunto de las variables independientes.

La colinealidad es un grave problema para el modelo multivariable, ya que las estimaciones obtenidas por el mismo son muy inestables, al hacerse más difícil separar el efecto de cada variable predictora.

Pues bien, para determinar si nuestro modelo sufre de colinealidad podemos construir una matriz donde se muestran los coeficientes de correlación, R, de unas variables con otras. En aquellos casos en los que observemos R altos, podremos sospechar que existe colinealidad. Ahora bien, si queremos cuantificar esto recurriremos a las otras dos patas del gato que hemos comentado al inicio: tolerancia y factor de inflación de la varianza.

Si elevamos el coeficiente R al cuadrado obtenemos el coeficiente de determinación (R²), que representa el porcentaje de la variación (o varianza) de una variable que es explicada por la variación en la otra variable. Así, nos encontramos con el concepto de tolerancia, que se calcula como el complementario de R² (1-R²) y que representa la proporción de la variabilidad de dicha variable que no se explica por el resto de las variables independientes incluidas en el modelo de regresión.

De esta forma, cuanto más baja sea la tolerancia, más probable será que exista colinealidad. Suele considerarse que existe colinealidad cuando R² es superior a 0,9 y, por tanto, la tolerancia está por debajo de 0,1.

Ya solo nos queda la tercera pata, que es el factor de inflación de la varianza. Este se calcula como el inverso de la tolerancia (1/T) y representa la proporción de la variabilidad (o varianza) de la variable que es explicada por el resto de las variables predictoras del modelo. Como es lógico, cuanto mayor sea el factor de inflación de la varianza, mayor será la probabilidad de que exista colinealidad. Generalmente se considera que existe colinealidad cuando el factor de inflación entre dos variables es mayor de 10 o cuando la media de todos los factores de inflación de todas las variables independientes es muy superior a uno.

Y aquí vamos a dejar los modelos multivariables por hoy. Ni que decir tiene que todo lo que hemos contado en la práctica se hace recurriendo a programas informáticos que nos calculan estos parámetros de manera sencilla.

Hemos visto aquí algunos de los aspectos de la regresión lineal múltiple, quizás el más utilizado de los modelos multivariables. Pero hay otros, como el análisis multivariante de la varianza (MANOVA), el análisis factorial o el análisis por conglomerados o clústeres. Pero esa es otra historia…

En busca de la causalidad

Publicado el 7 febrero, 2017 | por Manuel Molina | Deja un Comentario

En Medicina es frecuente que tratemos de buscar relaciones de causa efecto. Si queremos demostrar que el fármaco X produce un efecto, no tenemos más que tomar dos grupos de personas, a un grupo le damos el fármaco, al otro grupo no se lo damos y vemos si hay diferencias.

Pero la cosa no es tan sencilla, porque nunca podemos estar seguros de que las diferencias en efecto entre los dos grupos se deban en realidad a otros factores distintos al tratamiento que hemos empleado. Estos factores son los llamados factores de confusión, que pueden ser conocidos o desconocidos y que nos pueden sesgar los resultados de la comparación.

Para resolver este problema se inventó el elemento clave de un ensayo clínico, la aleatorización. Si repartimos los participantes en el ensayo entre las dos ramas de forma aleatoria conseguiremos que estas variables de confusión se repartan de forma homogénea entre las dos ramas del ensayo, con lo que cualquier diferencia entre las dos tendrá que ser debida a la intervención. Solo así podremos establecer relaciones de causa-efecto entre nuestra exposición o tratamiento y la variable de resultado que midamos.

El problema de los estudios cuasi-experimentales y de los observacionales es que carecen de aleatorización. Por este motivo, nunca podremos estar seguros de que las diferencias se deban a la exposición y no a cualquier variable confusora, por lo que no podemos establecer con seguridad relaciones causales.

Este es un inconveniente molesto, ya que muchas veces será imposible realizar ensayos aleatorizados ya sea por motivos éticos, económicos, de la naturaleza de la intervención o de lo que sea. Por eso se han inventado algunas argucias para poder establecer relaciones causales en ausencia de aleatorización. Una de estas técnicas es la de los propensity score que vimos en una entrada anterior. Otra es la que vamos a desarrollar hoy, que tiene el bonito nombre de regresión discontinua.

La regresión discontinua es un diseño cuasi-experimental que permite realizar inferencia causal en ausencia de aleatorización. Se puede aplicar cuando la exposición de interés se asigna, al menos parcialmente, según el valor de una variable aleatoria continua si esta variable cae por encima o por debajo de un determinado valor umbral. regresion-discontinua_umbral Pensemos, por ejemplo, en un fármaco hipocolesterolemiante que pautaremos cuando el colesterol LDL aumente por encima de un valor determinado, o de una terapia antirretroviral en un enfermo de sida que indicaremos cuando su contaje de CD4 disminuya por debajo de determinado valor. Existe una discontinuidad en el valor umbral de la variable que produce un cambio brusco en la probabilidad de asignación al grupo de intervención, tal como os muestro en la figura adjunta.

En estos casos en los que la asignación del tratamiento depende, al menos en parte, del valor de una variable continua, la asignación en las proximidades del umbral es casi como si fuese aleatoria. ¿Por qué? Porque las determinaciones están sujetas a una variabilidad aleatoria por error de muestreo (además de la propia variabilidad de las variables biológicas), lo que hace que los individuos que están muy cerca del umbral, por encima o por debajo, sean muy similares en cuanto a las variables que puedan actuar como confusoras (el estar por encima o por debajo del umbral puede depender de la variabilidad aleatoria del resultado de la medición de la variable), de manera similar a como ocurre en un ensayo clínico. A fin de cuentas, podemos pensar que un ensayo clínico no es más que un diseño de discontinuidad en el que el umbral es un número aleatorio.

La matemática de la regresión discontinua es solo para iniciados y no es mi intención explicarla aquí (primero tendría que entenderla yo), así que nos vamos a conformar con conocer algunos términos que nos servirán para entender los trabajos que empleen esta metodología.

La regresión discontinua puede ser nítida o difusa. En la nítida, la probabilidad de asignación cambia de cero a uno en el umbral (la asignación del tratamiento sigue una regla determinista). Por ejemplo, se inicia el tratamiento cuando se cruza el umbral, con independencia de otros factores. Por otra parte, en la difusa hay otros factores en juego que hacen que en el umbral la probabilidad de asignación cambie, pero no de cero a uno, sino que puede depender de esos otros factores añadidos.

Así, el resultado del modelo de regresión varía un poco según se trate de una regresión discontinua nítida o difusa. En el caso de la regresión nítida se calcula el llamado efecto causal medio, según el cual los participantes son asignados a la intervención con seguridad si traspasan el umbral. En el caso de la regresión difusa, la asignación ya no se realiza según un modelo determinista, sino según uno probabilístico (según el valor respecto al umbral y el de otros factores que el investigador puede considerar importantes). En estos casos hay que hacer un análisis por intención de tratamiento según la diferencia de la probabilidad de asignación cerca del punto de corte (algunos pueden no traspasar el umbral pero ser asignados a la intervención porque así lo considere el investigador según los otros factores).

Así, en el modelo probabilístico habrá que medir el efecto en los cumplidores (los asignados a la intervención), por lo que el modelo de regresión nos dará el efecto causal medio de los cumplidores, que es la medida típica de la regresión discontinua difusa.

Y creo que aquí lo vamos a dejar por hoy. No hemos hablado nada sobre la ecuación de regresión, pero baste decir que tiene en cuenta las pendientes de la función de probabilidad de asignación antes y después del umbral y una variable de interacción para la posibilidad de que los efectos del tratamiento sean heterogéneos a ambos lados del umbral. Como veis, todo bastante complicado, pero para eso están los paquetes estadísticos como R o Stata que implementan estos modelos sin apenas esfuerzo.

Para terminar, decir solo que lo habitual es ver modelos que utilizan regresión lineal para variables de resultado cuantitativas, pero existen extensiones del modelo que utilizan variables dicotómicas y técnicas de regresión logística, e incluso modelos con estudios de supervivencia y variables de tiempo a suceso. Pero esa es otra historia…

Censura

Publicado el 17 enero, 2017 | por Manuel Molina | Deja un Comentario

En el sentido más conocido de la palabra, censura es la acción de examinar una obra destinada al público, suprimiendo o modificando la parte que no se ajusta a determinados planteamientos políticos, morales o religiosos, para determinar si se puede o no publicar o exhibir. Entonces, ¿qué queremos decir en estadística cuando hablamos de datos censurados?. Nada que ver con la política, moral ni la religión. Para explicar lo que es un dato censurado tendremos que hablar primero de las variables de tiempo a suceso y de los análisis de supervivencia.

De manera general, podemos decir que hay tres tipos de variables: cuantitativas, cualitativas y de tiempo a suceso. Las dos primeras se entienden bastante bien en general, pero las de tiempo a suceso son un poco más complicadas de entender.

Imaginemos que queremos estudiar la mortalidad de esa terrible enfermedad que es la fildulastrosis. Podríamos contar el número de fallecidos al final del periodo del estudio y dividirlos por la población total al inicio. Por ejemplo, si al inicio hay 50 enfermos y se nos mueren cuatro durante el seguimiento, podríamos calcular la mortalidad como 4/50 = 0,08, o sea del 8%. Así, si hemos seguido a la población durante cinco años, podremos decir que la supervivencia de la enfermedad a los cinco años es del 92% (100-8 = 92).

Sencillo, ¿verdad? El problema es que esto solo es válido cuando todos los sujetos tienen el mismo periodo de seguimiento y no se producen pérdidas o abandonos a lo largo del estudio, situación que suele estar lejos de la realidad en la mayor parte de los casos.

En estos casos, lo correcto es medir no solo si se produce el fallecimiento (que sería una variable dicotómica), sino también cuándo se produce, teniendo en cuenta además el diferente periodo de seguimiento y las pérdidas. Así, utilizaríamos una variable de tiempo a suceso, que está compuesta por una variable dicotómica (el suceso que se mide) y una continua (el tiempo de seguimiento cuando se produce).

Siguiendo el ejemplo anterior, los participantes en el estudio podrían clasificarse en tres tipos: aquéllos que fallecen durante el seguimiento, los que permaneces vivos al final del estudio y los que se pierden durante el seguimiento.

De los que se mueren podemos calcular su supervivencia pero, ¿cuál es la supervivencia de los que están vivos al final del estudio? ¿Y cuál es la supervivencia de los que se pierden durante el seguimiento? Está claro que algunos de los perdidos pueden haber fallecido al final del estudio sin que nosotros lo detectemos, por lo que nuestra medida de la mortalidad no será exacta.

Y aquí es donde nos encontramos con los datos censurados. Todos aquellos que no presentan el evento durante un estudio de supervivencia se denominan censurados (las pérdidas y los que acaban el estudio sin presentar el evento). La importancia de estos datos censurados es que hay que tenerlos en cuenta al hacer el estudio de supervivencia, tal como veremos a continuación.

La metodología a seguir es confeccionar una tabla de supervivencia que tenga en cuenta los sucesos (en este caso las muertes) y los datos censurados, tal como vemos en la tabla adjunta.

Las columnas de la tabla representan lo siguiente: x, el número de año del seguimiento; Nx, el número de participantes vivos al inicio de ese año; Cx, el número de pérdidas de ese año (censurados); Mx, el número de fallecidos durante ese periodo; PM, probabilidad de morir en ese periodo; PSP, la probabilidad de sobrevivir en ese periodo (la probabilidad de no presentar el evento); y PSG, la probabilidad de supervivencia hasta ese momento. censura Como vemos, el primer año partimos de 50 participantes, de los cuales uno fallece. La probabilidad de fallecer en ese periodo es de 1/50 = 0,02, con lo que la probabilidad de supervivencia en el periodo (que es igual a la global por ser el primer periodo) es de 1-0,02 = 0,98.

En el segundo periodo partimos de 49 y no fallece ni se pierde nadie. La PM en el periodo es cero y la de supervivencia uno. Así, la probabilidad global será de 1×0,98 = 0,98.

En el tercer periodo seguimos con 49. Se pierden dos y fallece uno. La PM es de 1/49 = 0,0204 y la PSP de 1-0,0204 = 0,9796. Si multiplicamos la PSP por la global del periodo anterior, obtenemos la supervivencia global de este periodo: 0,9796×0,98 = 0,96.

En el cuarto periodo partimos de 46 participantes, produciéndose cinco pérdidas y dos fallecimientos. La PM será de 2/46 = 0,0434, la PSP de 1-0,0434 = 0,9566 y la PSG de 0,9566×0,96 = 0,9183.

Por último, en el quinto periodo partimos de 39 participantes. Tenemos dos censurados y ningún evento (fallecimiento). PM es cero, PSP es igual a uno (no se muere nadie en este periodo) y PSG 1×0,9183 = 0,9183.

Finalmente, teniendo en cuenta los datos censurados, podemos decir que la supervivencia global de la fildulastrosis es del 91,83% a los cinco años.

Y con esto vamos a dejarlo por hoy. Hemos visto cómo se construye una tabla de supervivencia con datos censurados para tener en cuenta el seguimiento desigual de los participantes y las pérdidas durante el seguimiento.

Solo dos reflexiones antes de terminar. En primer lugar, aunque se hable de análisis de supervivencia, el evento no tiene porqué ser el fallecimiento de los participantes. Puede ser cualquier evento que se produzca a lo largo del seguimiento del estudio.

En segundo lugar, las variables de tiempo a suceso y los datos censurados son la base para realizar otras técnicas estadísticas que estiman la probabilidad de producirse el evento en estudio en un momento determinado, como los modelos de regresión de Cox. Pero esa es otra historia…

Un caso de probabilidad engañosa

Publicado el 15 noviembre, 2016 | por Manuel Molina | Deja un Comentario

Hoy vamos a ver otro de esos ejemplos en los que la intuición sobre el valor de determinadas probabilidades nos juega malas pasadas. Y, para ello, vamos a utilizar nada menos que el teorema de Bayes, jugando un poco con las probabilidades condicionadas. Vamos a ver paso a paso cómo funciona.

¿Cuál es la probabilidad de que se produzcan dos sucesos? La probabilidad de que ocurra un suceso A es P(A) y la de que ocurra B, P(B). Pues bien, la probabilidad de que ocurran los dos es P(A∩B) que, si los dos sucesos son independientes, es igual a P(A) x P(B).

Imaginemos que tenemos un dado con seis caras. Si lo lanzamos una vez, la probabilidad de sacar, por ejemplo, un cinco es de 1/6 (un resultado entre los seis posibles). La de sacar un cuatro es, igualmente, 1/6. ¿Cuál será la probabilidad de sacar un cuatro, una vez que en la primera tirada sacamos un cinco?. Como las dos tiradas son independientes, la probabilidad de la combinación cinco seguida de cuatro será de 1/6 x 1/6 = 1/36.

Ahora pensemos otro ejemplo. Supongamos que en un grupo de 10 personas hay cuatro médicos, dos de los cuáles son cirujanos. Si tomamos uno al azar, la probabilidad de que sea médico es de 4/10 = 0,4 y la de que sea cirujano es de 2/10 = 0,2. Pero, si sacamos a uno y sabemos que es médico, la probabilidad de que sea cirujano ya no será de 0,2, porque los dos sucesos, ser médico y cirujano, no son independientes. Si es médico, la probabilidad de que sea cirujano será de 0,5 (la mitad de los médicos de nuestro grupo son cirujanos).

Cuando dos sucesos son dependientes, la probabilidad de que ocurran los dos será la probabilidad de ocurrir el primero, una vez que ocurre el segundo, por la probabilidad de ocurrir el segundo. Así que la P(médico∩cirujano) = P(cirujano|médico) x P(médico). Podemos generalizar la expresión de la siguiente manera:

P(A∩B) = P(A|B) x P(B), y cambiando de orden los componentes de la expresión, obtenemos la llamada regla de Bayes, de la siguiente forma:

P(A|B) = P(A∩B) / P(B).

La P(A∩B) será la probabilidad de B, una vez que se produce A, por la probabilidad de A = P(B|A) x P(A). Por otra parte, la probabilidad de B será igual a la suma de la probabilidad de producirse B una vez que se produzca A más la probabilidad de producirse B sin que ocurra A, lo que puesto de forma matemática queda de la siguiente forma:

P(B|A) x P(A) + P(B|A^c) x P(A^c), siendo P(A^c) la probabilidad de que no ocurra A.

Si sustituimos la regla inicial por sus valores desarrollados, obtendremos la expresión más conocida del teorema de Bayes:

$P(A|B)=\frac{P(B|A) \times P(A)}{P(B|A) \times P(A)+P(B|A^{{c}}) \times P(A^{{c}})}$ Vamos a ver cómo se aplica el teorema de Bayes con un ejemplo práctico. Pensemos en el caso de la fildulastrosis aguda, una grave enfermedad cuya prevalencia en la población es, afortunadamente, bastante baja, de uno por cada 1000 habitantes. Luego, la P(F) = 0,001.

Por suerte tenemos una buena prueba diagnóstica, con una sensibilidad del 98% y una especificidad del 95%. Supongamos ahora que yo me hago la prueba y me da un resultado positivo. ¿Tengo que asustarme mucho? ¿Cuál es la probabilidad de que realmente tenga la enfermedad? ¿Os parece que será alta o baja? Veámoslo.

Una sensibilidad del 98% quiere decir que la probabilidad de dar positivo cuando se tiene la enfermedad es de 0,98. Matemáticamente, P(POS|F) = 0,98. Por otra parte, una especificidad del 95% quiere decir que la probabilidad de que dé un resultado negativo estando sano es de 0,95. O sea, P(NEG|F^c) = 0,95. Pero nosotros lo que queremos saber no es ninguna de estas dos cosas, sino que realmente buscamos cuál es la probabilidad de estar enfermo una vez que damos positivo en la prueba, o sea, la P(F|POS).

Para calcularla, no tenemos más que aplicar el teorema de Bayes:

$P(F|POS)=\frac{P(POS|F) \times P(F)}{P(POS|F) \times P(F)+P(POS|F^{{c}}) \times P(F^{{c}})}$ A continuación, sustituimos los símbolos con sus valores y resolvemos la ecuación:

$P(F|POS)=\frac{0,98 \times 0,001}{0,98 \times 0,001+[(1-0,95) \times (1-0,001)]}=0,02$ Así que vemos que, en principio, no tengo que asustarme mucho cuando la prueba me da un resultado positivo, ya que la probabilidad de estar enfermo es solo de un 2%. Como veis, mucho más baja de lo que la intuición nos diría con una sensibilidad y una especificidad tan altas. ¿Por qué ocurre esto? Muy sencillo, porque la prevalencia de la enfermedad es muy baja. Vamos a repetir el experimento suponiendo ahora que la prevalencia es del 10% (0,1):

$P(F|POS)=\frac{0,98 \times 0,1}{0,98 \times 0,1+[(1-0,95) \times (1-0,1)]}=0,68$ Como veis, en este caso la probabilidad de estar enfermo si doy positivo sube hasta el 68%. Esta probabilidad es el conocido valor predictivo positivo que, como podemos comprobar, puede variar enormemente según la frecuencia del efecto que estemos estudiando.

Y aquí lo dejamos por hoy. Antes de terminar, dejadme advertiros que no busquéis qué es la fildulastrosis. Me sorprendería mucho que alguien la encontrase en algún libro de medicina. Además, tened cuidado de no confundir P(POS|F) con P(F|POS), ya que incurriríais en un pecado llamado falacia inversa o falacia de la transposición de los condicionales, que es un error grave.

Hemos visto como el cálculo de probabilidades se complica un poco cuando los sucesos no son independientes. También hemos aprendido lo poco de fiar que son los valores predictivos cuando cambia la prevalencia de la enfermedad. Por eso se inventaron los cocientes de probabilidades, que no dependen tanto de la prevalencia de la enfermedad que se diagnostica y permiten valorar mejor de forma global la potencia de la prueba diagnóstica. Pero esa es otra historia…

No te dejes llevar por los extremos

Publicado el 28 junio, 2016 | por Manuel Molina | Deja un Comentario

Ya vimos en una entrada anterior que los valores extremos de una distribución, los llamados outliers, pueden sesgar las estimaciones de los estadísticos que calculamos en nuestra muestra.

Un ejemplo típico es el de la media aritmética, que se desplaza en la dirección de los valores extremos, si los hay, tanto más cuanto más extremos sean los valores. Vimos que, para evitar este inconveniente, existían una serie de familiares de la media aritmética que se consideraban robustos o, lo que es lo mismo, que eran menos sensibles a la presencia de outliers. De todos estos, el más conocido es la mediana, aunque existen algunos más, como la media recortada, la winsorizada, la ponderada, la geométrica, etc.

Pues bien, algo parecido a lo que le pasa a la media ocurre también con la desviación típica, el estadístico de escala o dispersión utilizado con más frecuencia. La desviación típica o estándar también se ve sesgada por la presencia de valores extremos, obteniendo valores que son poco representativos de la dispersión real de la distribución.

Veamos el ejemplo que utilizábamos al hablar de los estimadores robustos de la media. Supongamos que medimos los valores de colesterol sérico en un grupo de personas y nos encontramos los siguientes valores (en mg/dl): 166, 143, 154, 168, 435, 159, 185, 155, 167, 152, 152, 168, 177, 171, 183, 426, 163, 170, 152 y 155. Como vemos, existen dos valores extremos (426 y 435 mg/dl) que nos sesgarán los estadísticos habituales que son la media y la desviación típica. En nuestro caso, podemos calcular la desviación típica y ver que su valor es de 83 mg/dl, claramente poco ajustado a la desviación de la mayoría de los valores respecto a cualquiera de las medidas de centralización robustas que podamos elegir.

¿Qué hacemos en este caso? Pues utilizar cualquiera de los estimadores robustos de la desviación, que hay varios. Algunos de ellos surgen a partir de los estimadores robustos de la media. Veamos algunos.

El primero, que surge a partir de la mediana, es la desviación absoluta mediana (DAM). Si recordáis, la desviación típica es la suma de las diferencias de cada valor con la media, elevadas al cuadrado, y dividida por el número de elementos, n (o por n-1 si lo que queremos es obtener un estimador no sesgado de la desviación típica poblacional). Pues bien, de modo similar, podemos calcular la mediana de las desviaciones absolutas de cada valor con la mediana de la muestra, según la siguiente fórmula

DAM = Mediana {|Xi – Me|}, para i=1 hasta n.

Podemos calcularla en nuestro ejemplo y vemos que vale 17,05 mg/dl, bastante más ajustado que la desviación típica clásica.

El segundo se calcula a partir de la media recortada. Esta, como su nombre indica, se calcula recortando un determinado porcentaje de la distribución, por sus extremos (la distribución tiene que estar ordenada de menor a mayor). Por ejemplo para calcular la media recortada al 20% de nuestro ejemplo quitaríamos un 10% por cada lado (dos elementos por cada lado: 143, 152, 426 y 435) y calcularíamos la media aritmética con los restantes. Pues bien, podemos calcular la desviación de la forma clásica con los elementos recortados, obteniendo el valor de 10,5 mg/dl.

Por último, en tercer lugar podríamos hacerlo siguiendo el razonamiento que se utiliza para calcular la media winsorizada. En este caso, en vez de eliminar los valores, los sustituiríamos por los valores más próximos sin eliminar. Una vez winsorizada la distribución, calculamos la desviación típica con los nuevos valores de la forma habitual. Su valor es de 9,3 mg/dl, similar a la anterior.

¿Cuál utilizamos de las tres?. Pues nos interesa utilizar una que se comporte de forma eficiente cuando la distribución sea normal (en estos casos la mejor es la desviación típica clásica) pero que no sea muy sensible cuando la distribución se aparte de la normal. En este sentido, la mejor es la desviación absoluta mediana, seguida de la desviación típica winsorizada muestral.

Un último consejo antes de finalizar. No os pongáis a calcular estas medidas a mano, ya que puede resultar muy laborioso, Los programas de estadística hacen los cálculos por nosotros sin el menor esfuerzo.

Y aquí terminamos. No hemos hablado nada de otros estimadores de la familia de los M-estimadores, como la varianza media biponderada o la varianza media de porcentaje ajustado. Estas medias son mucho más difíciles de comprender desde el punto de vista matemático, aunque son muy fáciles de calcular con el paquete informático adecuado. Pero esa es otra historia…

Una relación simple

Publicado el 31 mayo, 2016 | por Manuel Molina | Deja un Comentario

Hoy vamos a volver a hablar de la relación que puede existir entre dos variables. Vimos en una entrada anterior como podíamos medir la relación entre dos variables mediante el procedimiento de correlación, que nos medía la fuerza de relación entre dos variables cuando ninguna de las dos puede considerarse predictora de la otra. Esto es, cuando los valores de una no nos sirven para calcular los valores de la otra, aunque las dos varíen de una forma predecible.

Una cosa parecida, de la que vamos a hablar en esta entrada, es la regresión. Esta no solo explica la relación que hay entre dos variables, sino que podemos cuantificar cómo varía una de las variables, que llamaremos dependiente, con las variaciones de la otra variables, que será la independiente.

Pero todavía podemos llegar un paso más allá: los valores de la variable independiente nos pueden servir para predecir el correspondiente valor de la variable dependiente. Supongamos que medimos peso y talla y calculamos el modelo de regresión entre el peso y la talla. Si sabemos la talla de un individuo podemos utilizar la ecuación de regresión para estimar cuál será su peso (en este caso la talla es la variable independiente y el peso la dependiente).

Si llamamos x a la variable independiente e y a la variable dependiente, los modelos de regresión simple pueden representarse mediante la siguiente ecuación:

Función(y) = a + bx

En esta ecuación, a representa el valor de la función de y cuando x vale cero. Se suele llamar interceptor porque es el punto donde la representación gráfica de la recta de regresión cruza el eje de las y. Por su parte, b representa la llamada pendiente, que es la cantidad que varía y con las variaciones de x (si x aumenta en b unidades, y aumenta en b unidades).

¿Y qué significa función(y)?. Pues depende del tipo de variable que sea la variable dependiente. Sabemos que las variables se clasifican en cuantitativas (o continuas), cualitativas (nominales u ordinales) y de tiempo a suceso (también llamadas de supervivencia). Pues bien, según el tipo de la variable dependiente la función(y) será diferente porque aplicaremos un modelos de regresión simple diferente.

En el caso de variables continuas, el modelo de regresión que aplicamos es el de regresión lineal simple y la función de y será su media aritmética. La ecuación será la siguiente:

y = a + bx

Volviendo al ejemplo del peso y la talla, si sustituimos x por el valor de talla deseado y resolvemos la ecuación obtendremos el peso medio de los individuos de esa talla.

En el caso de que la variable dependiente sea cualitativa binaria utilizaremos un modelo de regresión logística. En este caso codificaremos la variable dependiente como cero y uno y la función de y ya no será la media, sino el logaritmo neperiano de la odds ratio del valor uno de la variable. Imaginemos que calculamos la relación entre peso (variable independiente) y sexo (variable dependiente). En este caso podríamos codificar como uno si es mujer y cero si es hombre, representando la recta de regresión de la siguiente forma:

Ln(OR) = a + bx

Si sustituimos x por el peso en cuestión y resolvemos la ecuación, obtendremos el logaritmo de la OR de ser mujer (el valor 1). Para obtener la OR debemos elevar el número e al resultado de la ecuación (hacer el antilogaritmo), obteniendo así la OR de que sea mujer. A partir de aquí es sencillo calcular el valor de la probabilidad de que sea mujer (p = OR/1+OR) u hombre (uno menos el valor de la probabilidad de que sea mujer).

Esta función del ln(OR) se expresa en muchas ocasiones como ln(p/1-p), ya que la odds ratio es la probabilidad de que un suceso ocurra (p) dividida de la probabilidad de que no ocurra (1-p). A esta función se la denomina logit, por lo que podemos ver también representada la regresión logística de la siguiente forma:

Logit(y) = a + bx

Por último, podemos encontrarnos el caso de que la variable dependiente sea una variable de tiempo a suceso. En este caso hay que utilizar un modelo de regresión de riesgos proporcionales de Cox. La estructura es muy similar a la de la regresión logística, solo que la función de y es el logaritmo de la hazard ratio en lugar del de la odds ratio:

Ln(HR) = a + bx

Igual que hacíamos con la regresión logística, para calcular el valor de la hazard ratio hay que hacer el antilogaritmo natural del producto de la ecuación de regresión (e elevado al resultado de la ecuación).

Y, aunque hay muchos más, estos son los tres modelos de regresión más utilizados. En todos estos casos hemos hablado de ecuaciones con una variable independiente, por lo que decimos que hablamos de regresión simple. Pero podemos meter todas las variables independientes que queramos, según la siguiente fórmula:

Función(y) = a + bx₁ + cx₂ + … + nx_n

Claro que ya no hablaríamos de regresión simple, sino de regresión múltiple, pero todo lo que hemos descrito sería igual de aplicable.

Y aquí lo vamos a dejar. Podríamos hablar del valor del interceptor y de la pendiente según la variable independiente sea continua o cualitativa, ya que se leen de forma un poco diferente. Pero esa es otra historia…

Ovejas negras

Publicado el 3 mayo, 2016 | por Manuel Molina | Deja un Comentario

Se dice que es una oveja negra aquél elemento de un grupo que va en dirección distinta o contraria a la del resto del grupo. Por ejemplo, en una familia de adictos a la telebasura, la oveja negra sería un miembro de esa familia que se desviviese por ver los documentales de la segunda cadena. Claro que si la familia es adicta a los documentales, la oveja negra se morirá por ver la telebasura. Siempre al revés.

En estadística hay algo parecido a las ovejas negras. Son los datos anómalos, también llamados datos extremos, pero más conocidos por su nombre en inglés: outliers.

Un outlier es una observación que parece inconsistente con el resto de los valores de la muestra, siempre teniendo en cuenta el modelo probabilístico supuesto que debe seguir la muestra. Como veis, es un dato que lleva la contraria a los demás, como una oveja negra.

El problema del outlier es que puede hacer mucho daño al estimar parámetros poblacionales a partir de una muestra. Vamos a recordar un ejemplo que vimos en otra entrada sobre el cálculo de medidas de centralidad robustas. Se trataba de un colegio con cinco maestros y un director fanático del futbol. Al hacer los contratos establece los siguientes sueldos: 1200 euros al mes para el profesor de ciencias, 1500 para el de mates, 800 para el de literatura y 1100 para el de historia. Pero resulta que se le antoja contratar a Zinedine Zidane como profesor de gimnasia, así que tiene que pagarle nada menos que 20000 euros mensuales.

¿Veis por dónde la va la cosa? Efectivamente, Zuzu es la oveja negra, el valor anómalo. Fijaos qué pasa si calculamos la media: 4920 euros al mes es el sueldo medio de los profesores de este centro. ¿Os parece una estimación real? Claramente no, el valor de la media está desplazada en la dirección del outlier, y se desplazaría más cuánto más extremo sea el valor anómalo. Si a Zuzu le pagasen 100000 euros, el sueldo medio ascendería a 20920 euros. Una locura.

Si un valor anómalo puede hacerle tanto daño a un estimador, imaginad lo que puede hacer con un contraste de hipótesis, en el que la respuesta es un aceptar o rechazar la hipótesis nula. Así que nos planteamos, ¿qué podemos hacer cuando descubrimos que entre nuestros datos hay una (o varias) ovejas negras? Pues podemos hacer varias cosas.

La primera que se nos pasa por la cabeza es tirar el outlier a la basura. Prescindir de él a la hora de analizar los datos. Esto estaría bien si el valor extremo es producto de un error en la recogida de los datos pero, claro, podemos prescindir de datos que dan información adicional. En nuestro ejemplo, el outlier no es ningún error, sino que es producto del historial deportivo del profesor en cuestión. Necesitaríamos algún método más objetivo para poder decidir suprimir el outlier, y aunque existen unas pruebas llamadas de discordancia, tienen sus problemas.

La segunda cosa que podemos hacer es identificarlo. Esto significa que tenemos que averiguar si el valor es tan extremo por alguna razón concreta, como pasa en nuestro ejemplo. Un valor extremo puede estar señalando algún hallazgo importante y no tenemos porqué desdeñarlo con rapidez, sino tratar de interpretar su significado.

En tercer lugar, podemos incorporarlos. Como hemos dicho al definirlos, el outlier lleva la contraria a los demás datos de la muestra según el modelo de probabilidad que suponemos que sigue la muestra. A veces, un dato extremo deja de serlo si asumimos que los datos siguen otro modelo. Por ejemplo, un outlier puede serlo si consideramos que los datos siguen una distribución normal pero no si consideramos que siguen una logarítmica.

Y, en cuarto lugar, la opción más correcta de todas: utilizar técnicas robustas para hacer nuestras estimaciones y nuestros contrastes de hipótesis. Se llaman técnicas robustas porque se afectan menos por la presencia de valores extremos. En nuestro ejemplo con los profesores utilizaríamos una medida de centralidad robusta como es la mediana. En nuestro caso es de 1200 euros, bastante más ajustada a la realidad que la media. Además, aunque le paguen a Zuzu 100000 euros al mes, la mediana seguirá siendo de 1200 euros mensuales.

Y con esto terminamos con los valores anómalos, esas ovejas negras que se mezclan con nuestros datos. No hemos comentado nada por simplificar, pero también podríamos tratar de averiguar cómo afecta el outlier a la estimación del parámetro, para lo cual existe toda una serie de metodología estadística basada en la determinación de la llamada función de influencia. Pero esa es otra historia…

No abuses de las tartas

Publicado el 15 marzo, 2016 | por Manuel Molina | Deja un Comentario

¡Qué ricas las tartas! El problema es que, como ya sabéis, lo que no está mal visto socialmente, o engorda o produce cáncer. Y las tartas no podían ser menos, así que hay que procurar no comer demasiado para que no se nos vayan al michelín o a otros sitios peores.

Pero hay una tarta que no engorda nada en absoluto (tampoco produce cáncer) y es el diagrama de tarta, que se utiliza con mucha frecuencia en estadística. ¿He dicho con mucha frecuencia? Probablemente me quede corto. Como no engorda ni tiene otros efectos perjudiciales para la salud hay tendencia a abusar de su uso.

tarta El gráfico de tarta, cuyo nombre correcto es gráfico de sectores, es muy sencillo de dibujar. Consiste en una circunferencia cuya área representa el total de los datos. Así, a cada categoría se le asigna un área que será directamente proporcional a su frecuencia. De esta forma, las categorías más frecuentes tendrán áreas mayores, de modo que de un vistazo podemos hacernos una idea de cómo se distribuyen las frecuencias en las categorías.

Hay tres formas de calcular el área de cada sector. La más sencilla es multiplicar la frecuencia relativa de cada categoría por 360°, obteniendo los grados de ese sector.

La segunda es utilizando la frecuencia absoluta de la categoría, según la siguiente regla de tres:

$\frac{Frecuencia\ absoluta}{Frecuencia\ total\ de\ datos}=\frac{Grados\ del\ sector}{360}$

Por último, la tercera forma consiste en utilizar las proporciones o porcentajes de las categorías:

$\frac{%\ de\ la\ variable}{100%}=\frac{Grados\ del\ sector}{360}$

Las fórmulas son muy sencillas pero, de todas formas, no habrá necesidad de recurrir a ellas porque el programa con el que hagamos el gráfico lo hará por nosotros.

El gráfico de sectores está diseñado para representar variables categóricas nominales, aunque no es raro ver tartas representando variables de otros tipos. Sin embargo, y en mi humilde opinión, esto no es totalmente correcto.

Por ejemplo, si hacemos un gráfico de sectores para una variable cualitativa ordinal estaremos perdiendo la información sobre la jerarquía de las variables, por lo que resultaría más correcto utilizar un gráfico que permita ordenar las categorías de menos a más. Y este gráfico no es otro que el diagrama de barras.

El diagrama de sectores será especialmente útil cuando haya pocas variables. Si hay muchas la interpretación deja de ser tan intuitiva, aunque siempre podemos completar el gráfico con una tabla de frecuencias que nos ayude a interpretar mejor los datos. Otro consejo es tener mucho cuidado con los efectos en 3D a la hora de dibujar las tartas. Si nos pasamos de elaborados el gráfico perderá claridad y será más difícil de leer.

Para terminar, deciros que tampoco tiene sentido utilizar una tarta para representar una variable cuantitativa. Para eso existe otro procedimiento más adecuado, que es el de utilizar un histograma, gráfico que mejor representa la distribución de frecuencias de una variable cuantitativa continua. Pero esa es otra historia…

¿Por qué sobra uno?

Publicado el 2 febrero, 2016 | por Manuel Molina | Deja un Comentario

Hoy vamos a hablar sobre uno de esos misterios de la estadística que muchos desconocen por qué son cómo son. Me refiero a si dividir entre n (el tamaño muestral) o entre n-1 para calcular las medidas de centralización y dispersión de una muestra, concretamente su media (m) y su desviación estándar (s).

La media sabemos todos lo que es. Su propio nombre lo dice, es el promedio de valores de una distribución de datos. Para calcularla sumamos todos los valores de la distribución y dividimos entre el total de elementos, o sea, entre n. Aquí no hay duda, dividimos entre n y obtenemos la medida de centralización más utilizada.

Por su parte, la desviación estándar, es una medida de la desviación media de cada valor respecto a la media de la distribución. Para obtenerla calculamos las diferencias de cada elemento con la media, las elevamos al cuadrado para que las negativas no se anulen con las positivas, las sumamos, las dividimos entre n y, por último, obtenemos la raíz cuadrada. Al ser la media de cada desviación, habrá que dividir las sumas de las desviaciones entre el total de elementos, n, como hacíamos con la media, según la conocida fórmula de la desviación estándar.

Sin embargo, en muchas ocasiones vemos que, para calcular la desviación estándar, dividimos entre n-1. ¿Por qué nos sobra un elemento?. Veámoslo.

Nosotros habitualmente trabajamos con muestras, de las que obtenemos sus medidas de centralización y dispersión. Sin embargo, lo que a nosotros nos interesaría saber en realidad es el valor de los parámetros en la población de la que procede la muestra. Por desgracia, no podemos calcular estos parámetros directamente, pero sí que podemos estimarlos a partir de los estadísticos de la muestra. Así, queremos saber si la media de la muestra, m, es un buen estimador de la media de la población, µ. Además, queremos saber si la desviación estándar de la muestra, s, es un buen estimador de la desviación de la población, que llamaremos σ.

Vamos a hacer un experimento para ver si m y s son buenos estimadores de µ y σ. Para ello vamos a utilizar el programa R. Os dejo el listado de comandos (script) en la figura adjunta por si queréis reproducirlo conmigo.

Primero generamos una población de 1000 individuos con una distribución normal con media de 50 y desviación estándar de 15 (µ = 50 y σ = 15). Una vez hecho, vamos a ver primero qué pasa con la media.

Si obtenemos una muestra de 25 elementos de la población y calculamos su media, esta se parecerá a la de la población (siempre que la muestra sea representativa de la población), pero puede haber diferencia debidas al azar. Para soslayar estas diferencias, obtenemos 50 muestras diferentes, con sus 50 medias. Estas medias siguen una distribución normal (la llamada distribución de muestreo), cuya media es la media de todas las que hemos obtenido de las muestras. Si extraemos 50 muestras con R y hallamos la media de sus medias, vemos que esta vale 49,6, lo que es casi igual a 50. Vemos, pues, que con las medias de las muestras podemos estimar bien el valor de la media de la distribución.

¿Y qué pasa con la desviación estándar? Pues si hacemos lo mismo (extraer 50 muestras, calcular su s y, por último, calcular la media de la 50 s) obtenemos una s media de 14,8. Esta s es bastante próxima al valor 15 de la población, pero se ajusta menos que el valor de la media. ¿Por qué?

La respuesta es que la media muestral es lo que se llama un estimador no sesgado de la media poblacional, ya que el valor medio de la distribución de muestreo es un buen estimador del parámetro en la población. Sin embargo, con la desviación estándar no pasa lo mismo, porque es un estimador sesgado. Esto es así porque la variación de los datos (que es a fin de cuentas lo que mide la desviación estándar) será mayor en la población que en la muestra, al tener la población un tamaño mayor (a mayor tamaño, mayor posibilidad de variación). Por eso dividimos por n-1, para que el resultado sea un poco más alto.

Si hacemos el experimento con R dividiendo entre n-1 obtenemos una desviación estándar no sesgada de 15,1, algo más próxima que la que obteníamos dividiendo entre n. Este estimador (dividiendo entre n-1) sería un estimador no sesgado de la desviación estándar poblacional. Entonces, ¿cuál empleamos? Si queremos saber la desviación estándar de la muestra podemos dividir entre n, pero si lo que queremos es una idea de cuánto vale el valor teórico en la población, el estimador se aproximará más al valor de σ si dividimos entre n-1.

Y aquí terminamos este galimatías. Podríamos hablar de cómo podemos obtener no solo el estimador a partir de la distribución de muestreo, sino también su intervalo de confianza, que nos diría entre que valores está el parámetro de la población, con un nivel de confianza determinado. Pero esa es otra historia…

Una caja con bigotes

Publicado el 24 noviembre, 2015 | por Manuel Molina | Deja un Comentario

No me negaréis que es un nombre bastante curioso para un gráfico. Al pensar en el nombre a uno se le viene a la cabeza una imagen de una caja de zapatos con bigotes de gato, pero, en realidad, yo creo que este gráfico se parece más a las naves de caza de tipo tie fighter de la Guerra de las Galaxias.

En cualquier caso, el gráfico de caja y bigotes, cuyo nombre formal es el de gráfico de caja (boxplot en inglés), es empleado con muchísima frecuencia en estadística por sus interesantes capacidades descriptivas.

Para saber de qué hablamos, tenéis representados dos gráficos de caja en la primera figura que os adjunto. Como veis, el gráfico, que puede representarse en vertical y en horizontal, consta de una caja y dos segmentos (los bigotes).

Describiendo la representación vertical, que quizás sea la más habitual, el borde inferior de la caja representa el percentil 25º de la distribución o, lo que es lo mismo, el primer cuartil. Por su parte, el borde superior (que se corresponde con el borde derecho de la representación horizontal) representa el percentil 75º de la distribución o, lo que es lo mismo, el tercer cuartil. De esta manera, la amplitud de la caja se corresponde con la distancia entre los percentiles 25º y 75º, que no es otra que el recorrido o rango intercuartílico. Por último, en el interior de la caja hay una línea que representa la mediana (o segundo cuartil) de la distribución. A veces puede verse una segunda línea que representa la media, aunque no es lo más habitual.

Vamos ahora con los bigotes. El superior se extiende hasta el valor máximo de la distribución, pero no puede llegar más allá de 1,5 veces el rango intercuartílico. Si existen valores más altos que la mediana más 1,5 veces el rango intercuartílico, éstos se representan como puntos más allá del extremo del bigote superior. Estos puntos son los denominados valores extremos, outliers en inglés. Vemos en nuestro ejemplo que hay un outlier que se sitúa más allá del bigote superior. Si no hay valores extremos u outliers, el máximo de la distribución lo marca el extremo del bigote superior. Si los hay, el máximo será el valor extremo más alejado de la caja.

Por añadidura, todo esto vale para el bigote inferior, que se extiende hasta el valor mínimo cuando no hay valores extremos o hasta la mediana menos 1,5 veces el rango intercuartílico cuando los haya. En estos casos, el valor mínimo será el outlier más alejado de la caja por debajo del bigote inferior.

Pues ya podemos comprender la utilidad del gráfico de caja. De un vistazo podemos obtener la mediana y el rango intercuartílico de su distribución e intuir la simetría de la distribución. Es fácil imaginarse cómo es el histograma de una distribución viendo su gráfico de caja, como podéis ver en la segunda figura. El primer gráfico corresponde a una distribución simétrica, próxima a la normal, ya que la mediana está centrada en la caja y los dos bigotes son más o menos simétricos.Si nos fijamos en la distribución central, la mediana está desplazada hacia el borde inferior de la caja y el bigote superior es más largo que el inferior. Esto es así porque la distribución tiene la mayoría de los datos hacia la izquierda y una larga cola hacia la derecha, como puede verse en su histograma. Lo mismo que hemos dicho para la distribución central vale parta la tercera, pero en este caso el bigote largo es el inferior y el sesgo es hacia la izquierda.

Por último, este tipo de gráfico sirve también para comparar varias distribuciones. En la tercera imagen que os adjunto podéis ver dos distribuciones aparentemente normales y con medianas muy similares. Si queremos hacer un contraste de hipótesis sobre la igualdad de sus medias, primero tenemos que saber si sus varianzas son iguales (si existe homocedasticidad) para saber qué tipo de test hay que utilizar.

Si comparamos las dos distribuciones, vemos que la amplitud de la caja y de los bigotes es mucho mayor en la primera que en la segunda, por lo que podemos concluir que la varianza de la primera distribución es mucho mayor, por lo que no podremos asumir la igualdad de varianzas y tendremos que aplicar la corrección pertinente.

Y esto es todo lo que quería contar sobre esta caja con bigotes, que tan útil resulta en estadística descriptiva. Ni que decir tiene que, aunque nos sirve para saber aproximadamente si la distribución se ajusta a una normal o si las varianzas de varias distribuciones son semejantes, existen pruebas específicas para estudiar estos puntos de forma matemática. Pero esa es otra historia…