Estadística con múltiples variables e introducción al lenguaje R

Algunos Conceptos Básicos de Estadística

Introducción

En este capítulo se presentarán algunos conceptos de estadística básica de una variable que serán necesarios para el desarrollo de métodos en varias variables. Se presentarán algunos conceptos de inferencia estadística, junto con algunas distribuciones elementales. Suponemos que el lector ha tenido algún contacto con los conceptos de estadística en una variable como parte de algún curso a nivel universitario.

Variables Aleatorias

Los modelos matemáticos se dividen en tres clases generales: (1) Puramente deterministas, (2) estáticos o deterministas con algunos componentes aleatorios sencillos, y (3) estocásticos. Cualquier observación de un modelo determinista es estrictamente una función de sus parámetros y variables tales como el tiempo, el espacio, o energías o estímulos varios. Los modelos puramente deterministas son típicos de la física clásica, donde por ejemplo, la distancia recorrida por un cuerpo en caída libre es una función de su altura y velocidad inicial. Si se ignoran efectos transitorios como la turbulencia atmosférica, error observacional entre otros, el desplazamiento se puede calcular exactamente dados un intervalo de tiempo y la constante gravitacional.

En el segundo tipo de modelo, cada observación es una función de un componente estrictamente determinista y un término aleatorio atribuible a errores en las medidas o a las variaciones en el muestreo de las observaciones o de las variables de entrada. Los modelos examinados aquí son principalmente de este tipo con la restricción de que la componente aleatoria será simplemente agregada a la parte determinista.

Los modelos estocásticos se construyen a partir de eventos o componentes fundamentalmente aleatorios para explicar fenómenos dinámicos o evolutivos: pueden ser tan simples como una secuencia de lanzamientos de una moneda al aire o tan complejos como los procesos de nacimiento y muerte que describen la dinámica de una población biológica.

Definamos ahora con más precisión que significa el concepto de variación aleatoria o los componentes aleatorios en los modelos tipo 2 y 3 mencionados anteriormente. Empezamos definiendo una variable aleatoria discreta o una que pueda contener un número contable de valores.

Supóngase que un experimento sólo puede resultar en uno y sólo uno de k resultados llamados $E_{1},...,E_{k}$. Estos resultados son mutuamente exclusivos, en el sentido de que la ocurrencia de uno niega la existencia de cualquier otro. A cada evento se le asigna un valor numérico $p_{i}$ entre cero y uno. Este valor numérico se le conoce como la probabilidad $P(E_{i})$ de ese evento. Entonces $p_{i}$ es la probabilidad de que ocurra el evento $E_{i}$ en una corrida del experimento. Dentro del marco del experimento se le asigna la probabilidad de cero a los eventos imposibles y la probabilidad de uno a los eventos que ocurrirán con entera certidumbre. De modo que, por exclusividad mutua de los eventos, en una corrida del experimento

${\begin{matrix}P(E_{i}\cap E_{j})&=&0\\P(E_{i}\cup E_{j})&=&p_{i}+p_{j}\end{matrix}}$ donde el símbolo de intersección $\cap$ denota el evento $E_{i}$

y $E_{j}$ y el símbolo de unión $\cup$ indica el evento $E_{i}$ o $E_{j}$.

Usando la propiedad aditiva de las probabilidades de eventos mutuamente exclusivos la probabilidad total del conjunto de eventos es

$P(E_{1}\cup E_{2}\cdots \cup E_{k})=p_{1}+\cdots +p_{k}=1$

Ahora se le asigna el valor numérico $x_{i}$ al i-ésimo resultado, donde, por conveniencia, los resultados se ordenan ascendentemente según sus valores $x_{i}$. Se define entonces la variable aleatoria discreta $X$ como la cantidad que toma el valor $x_{i}$ con probabilidad de $p_{i}$ en cada corrida del experimento. Por ejemplo, si el experimento consiste en lanzamientos de una moneda, el resultado "cara" puede tener una cierta probabilidad $p_{1}$ y el resultado "cruz" puede tener una probabilidad $p_{2}$. Esta variable aleatoria puede ser descrita entonces por una función de probabilidad $f(x_{i})$ la cual especifica las probabilidades con las que la variable aleatoria $X$ adquiere los valores de cara o cruz.

$X$	$f(x_{i})$
cara	$p_{1}$
cruz	$p_{2}$

Se hace notar que siempre la suma de todas las probabilidades debe ser igual a uno y que implícitamente se le ha asignado una probabilidad de cero a eventos irrelevantes tales como la moneda cayendo al canto.

Las variables aleatorias que se verán más adelante tomarán valores sobre alguna región continua en lugar de hacerlo sobre un conjunto de eventos discretos y se les llamará variables aleatorias continuas, variantes continuas o simplemente variantes. Todos estos términos son sinónimos. La variable aleatoria continua $X$ definida sobre el dominio de los números reales se caracteriza por su función de distribución

$F(x)=P(X\leq x)\;\;-\infty <x<\infty$

la cual proporciona la probabilidad de que $X$ sea menor o igual que algun valor de $x$ en su dominio. Ya que $X$ es continua, $P(X=x)=0$.

Si $F(x)$ es una función absolutamente continua, la versión continua de la función discreta de probabilidad es la función de densidad

$f(x)={\dfrac {dF(x)}{dx}}$

Al mismo tiempo, por la propiedad de continuidad absoluta

$F(x)=\int _{-\infty }^{x}f(u)\;du$

y de esta definición se genera el término equivalente función acumulativa de distribución de $F(x)$. Nótese que éstas definiciones son perfectamente generales: si la variable aleatoria se encuentra definida en algún intervalo de la recta de los números reales, fuera de ese intervalo $f(x)$ es igual a cero, y $F(x)$ es igual a cero y a uno a la izquierda y a la derecha del intervalo, respectivamente.

Cuando se consideran con respecto a su función de densidad $f(x)$, los valores en el intervalo se definen como la población o el universo de la variable aleatoria $X$.

Las propiedades de una variable aleatoria se visualizan comúnmente en términos de su función de densidad. Los términos exponencial o rectangular se refieren a estas representaciones. La figura muestra las funciones de densidad de tres variables conocidas. Las funciones de distribución de las variables exponenciales y normales se obtienen integrando directamente. La función de distribución normal sólo se puede evaluar por integración numérica.

Las funciones de la figura involucran parametros que determinan sus posiciones y sus formas. Los parámetros ${\frac {1}{2}}(a+b)$ de la densidad rectangular y $\mu$ de la función normal son parámetros de ubicación porque sus valores especifican las posiciones de las densidades en el eje real. El intervalo $b-a$ de la variable rectangular, el parámetro exponencial $\beta$, y el parámetro $\sigma$ de la densidad normal, son parámetros de escala porque cambios en sus valores equivalen a cambios en las unidades de las variables. Los incrementos en los valores de estos parámetros implican un mayor extendimiento de la función de densidad y por lo tanto más variación en la variable aleatoria. En general, si la densidad $f(x;\alpha ,\beta )$ se puede expresar como

$f\left({\dfrac {x-\alpha }{\beta }}\right)$

entonces $\alpha$ es un parámetro de ubicación y $\beta$ es un parámetro de escala. Se hace notar que todas las funciones de densidad de esta forma contienen el factor $1/\beta$ asociado con el elemento diferencial $dx$.

Las variables aleatorias y sus densidades tienen una interpretación física: considérese la función de densidad $f(x)$ de la variable aleatoria $X$ como la función que define la densidad de una barra sólida que se encuentra posicionada en el eje de las $x$. El k-ésimo momento alrededor del origen es

$\mu _{k}^{'}=\int _{-\infty }^{\infty }x^{k}f(x)\;dx$

El primer momento es la media o el valor esperado

$E(X)=\int _{-\infty }^{\infty }xf(x)\;dx$

de la variable aleatoria y corresponde al concepto físico del centro de gravedad de la vara. El símbolo $E$ denota la operación de calcular el valor esperado, y se conoce como el operador de esperanza. Si $c$ y $k$ son cantidades no aleatorias, entonces podemos aseverar las siguientes propiedades de las esperanzas:

${\begin{aligned}E(c)&=c\\E(cX)&=cE(X)\\E(k+cX)&=k+cE(X)\end{aligned}}$

Como ejemplo, considérese la función de densidad rectangular de la figura

Figura 1 - Funciones de Densidad Rectangular y Exponencial

. Por definición:

${\begin{aligned}E(X)&={\dfrac {1}{b-a}}\int _{a}^{b}x\;dx\\&={\dfrac {1}{b-a}}{\dfrac {1}{2}}x^{2}{\Bigr |}_{a}^{b}\\&={\frac {1}{2}}(b+a)\end{aligned}}$

Se puede intuir que el valor medio es el punto medio de los límites de la función de densidad, y si se desea, se puede reparametrizar la función de densidad en términos de la media como parámetro de ubicación y el intervalo $w=b-a$ como el parámetro de escala. Entonces, la nueva densidad será:

$h(x)={\begin{cases}0,&-\infty <x<E(X)-{\frac {1}{2}}w\\{\frac {1}{w}},&E(X)-{\frac {1}{2}}w\leq x\leq E(X)+{\frac {1}{2}}w\\0,&E(X)+{\frac {1}{2}}w<x<\infty \end{cases}}$

La varianza de una variable aleatoria es el valor esperado de las desviaciones cuadradas alrededor de la media, conocido también como el segundo momento central. Se denota la varianza de $X$ como

${\text{var}}(X)$

o también como $\sigma _{x}^{2}$

Por definición

${\begin{aligned}{\text{var}}(X)=&\int _{-\infty }^{\infty }\left[x-E(X)\right]^{2}f(x)\;dx\\=&\int _{-\infty }^{\infty }x^{2}f(x)\;dx-\left[E(X)\right]^{2}\\=&E(X^{2})-\left[E(X)\right]^{2}\end{aligned}}$

y el simbolo var aplicado a cualquier variable aleatoria denotará el cálculo de la varianza. Si $c$ y $k$ son constantes, la varianza posee las siguientes propiedades:

${\begin{aligned}{\text{var}}(c)&=0\\{\text{var}}(cX)&=c^{2}{\text{var}}(X)\\{\text{var}}(X+c)&={\text{var}}(X)\end{aligned}}$

La primera propiedad simplemente afirma que una variable no aleatoria tiene cero varianza. La naturaleza cuadrática de la varianza se encuentra reflejada en la segunda propiedad ya que un cambio de escala $c$ altera la varianza en $c^{2}$. El tercer resultado proviene de la definición y afirma que la varianza no se ve afectada al cambiar el origen en el eje de las $X$. Se puede demostrar que las varianzas de las funciones de densidad rectangular, exponencial, y normal

Figura 2 - Funciones de Densidad Rectangular, Exponencial, y Gaussiana

son ${\frac {1}{12}}(b-a)^{2}$, $\beta ^{2}$, y $\sigma ^{2}$ respectivamente.

En muchas ocasiones es deseable tener una medida de dispersión que posea las mismas unidades de la variable. La desviación estándar de $X$ es la raiz cuadrada positiva de la varianza:

$\sigma _{x}=+{\sqrt {{\text{var}}(X)}}$

Hacemos notar que los parámetros naturales de la densidad normal en la Figura 2 son la media y la desviación estándar.

Variables Independientes

En las primeras partes de este capítulo nos referimos vagamente a las variables aleatorias "independientes". Ahora proporcionaremos una definición precisa de la independencia de variables. Como se verá en el Capitulo 3, se puede extender el concepto de distribución y de funciones de densidad a múltiples variables, y si se puede expresar la función de distribución conjunta de $X_{1},\ldots ,X_{p}$ como

$F(x_{1},\ldots ,x_{p})=\int _{-\infty }^{x_{p}}\cdots \int _{-\infty }^{x_{1}}f(u_{1},\ldots ,u_{p})du_{1}\cdots du_{p}$

donde $f(u_{1},\ldots ,u_{p})$ es la densidad conjunta, entonces las variables se dice que son independientes si y solo si

$F(x_{1},\ldots ,x_{p})=F_{1}(x_{1})\cdots F_{p}(x_{p})$

donde $F_{i}(x_{i})$ es la función de distribución de la variable sencilla $X_{i}$. Alternativamente, la independencia es válida si y sólo si la factorización de

$f(x_{1},\ldots ,x_{p})=f_{1}(x_{1})\cdots f_{p}(x_{p})$

de la densidad conjunta en el producto de las funciones individuales de densidad es válida.

El momento de productos

$E(X_{1}^{k_{1}}\cdots X_{p}^{k_{p}})=\int _{-\infty }^{\infty }\cdots \int _{-\infty }^{\infty }x_{1}^{k_{1}}\cdots x_{p}^{k_{p}}f(x_{1},\ldots ,x_{p})\;dx_{1}\cdots dx_{p}$

de las variables $X_{1},\ldots ,X_{p}$ se factoriza como el producto

$E(X_{1}^{k_{1}})\cdots E(X_{p}^{k_{p}})$

de los $k$ momentos individuales si las variables son independientes. De este resultado se puede concluír que la varianza de la suma de variables aleatorias independientes es simplemente la suma de las varianzas individuales: ${\text{var}}(X_{1}+\cdots +X_{p})={\text{var}}(X_{1})+\cdots +{\text{var}}(X_{p}).$

Variables Aleatorias Normales (Gaussianas)

En esta sección se describirán algunas propiedades de una sola variable aleatoria normal como preparación para los resultados y técnicas que se basarán en la distribución normal multivariable. Recuérdese que la función de densidad normal es

$f(x)={\dfrac {1}{\sqrt {2\pi \sigma }}}\exp \left[-{\dfrac {1}{2}}\left({\dfrac {x-\mu }{\sigma }}\right)^{2}\right]\;\;-\infty <x<\infty$

y que su función de distribución se encuentra dada por la integral

$P(x)={\dfrac {1}{\sqrt {2\pi \sigma }}}\int _{-\infty }^{x}\exp \left[-{\dfrac {1}{2}}\left({\dfrac {x-\mu }{\sigma }}\right)^{2}\right]dx={\dfrac {1}{\sqrt {2\pi }}}\int _{-\infty }^{(x-\mu )/\sigma }\exp(-{\frac {1}{2}}u^{2})du=\Phi \left({\dfrac {x-\mu }{\sigma }}\right)$

donde $\Phi (z)$ denota la función de distribución normal estándar con media igual a cero y varianza igual a uno. Se denotará la distribución de una variable normal por medio del llamado símbolo de Wilks: $N(\mu ,\sigma ^{2})$; la distribución estandar será indicada por $N(0,1)$.

Los valores de $\Phi (z)$ se pueden encontrar fácilmente en cualquier tabla de funciones estadísticas. Por otro lado, el porcentaje 100α de la distribución normal unitaria se define como el valor $z_{\alpha }$ tal que

$\alpha =P(Z>z_{\alpha })=1-\Phi (z_{\alpha })$

Estas relaciones entre porcentajes y sus probabilidades se ilustran en la figura 3

Fig 3-Relaciones entre funciones de densidad y probabilidades

para la densidad normal y la función de distribución.

Ahora, permítasenos ofrecer la justificación donde mucha de la metodología estadística se basa para justificar la suposición de una población normal. Existe el llamado teorema del límite central el cual afirma que las variables aleatorias que son la suma de muchos efectos independientes tienden a encontrarse normalmente distribuídas a medida que el número de efectos va creciendo. Formalmente:

Si las variables aleatorias $X_{1},\ldots ,X_{N}$ se encuentran distribuídas independientemente de acuerdo a alguna función común de distribución con media $E(X_{i})=\mu$ y varianza finita ${\text{var}}(X_{i})=\sigma ^{2}$, entonces a medida que el número de variantes $N$ aumenta sin cotas, la variante

$Z={\dfrac {\sum _{i=1}^{N}X_{i}-N_{\mu }}{\sqrt {N\sigma ^{2}}}}$

converge en su distribución a una variante normal con media cero y varianza igual a uno.

Se puede encontrar una demostración de ésta versión del teorema del límite central en cualquier texto de probabilidad o estadística matemática. Una consecuencia inmediata del teorema es que la función ${\sqrt {N}}({\bar {X}}-\mu )$ de la media de las muestras de una secuencia de variables aleatorias independientes

${\bar {X}}={\frac {1}{N}}(X_{i}+\ldots +X_{N})$

cuya distribución común tiene media $\mu$ y varianza $\sigma ^{2}$ tiende a distribuírse como una variante de tipo $N(0,\sigma ^{2})$ a medida que $N$ se vuelve más grande. Una y otra vez estaremos invocando una propiedad muy importante de la distribución normal: Si $X_{1},\ldots ,X_{n}$ son variantes independientes con distribuciones $N(\mu _{1},\sigma _{1}^{2}),\ldots ,N(\mu _{n},\sigma _{n}^{2})$ entonces la combinación lineal

$Y=a_{1}X_{1}+\ldots +a_{n}X_{n}$

también se encuentra normalmente distribuída con media $\sum _{i=1}^{n}a_{i}\mu _{i}$ y varianza $\sum _{i=1}^{n}a_{i}^{2}\sigma _{i}^{2}$. Veremos más adelante que este resultado se puede extender a una combinación lineal de variantes normales dependientes.

Nota: Muchas de las siguientes funciones hacen uso de la función gamma $\left(\Gamma (x)\right)$. En el recuadro al final de esta sección se hace una breve presentación de esta función y sus propiedades.

La distribución ji cuadrada

Muchas distribuciones se pueden derivar a partir de diferentes transformaciones de un conjunto de variantes normalmente distribuídas. Una de las más importantes es la de la variante ji cuadrada. Si las variantes $X_{1},\ldots ,X_{n}$ son independientes y se encuentran normalmente distribuídas con media cero y varianza igual a uno, entonces

$\chi ^{2}=X_{1}^{2}+\ldots +X_{n}^{2}$

tiene la función de densidad

$f(\chi ^{2})={\frac {2^{-n/2}}{\Gamma (n/2)}}(\chi ^{2})^{(n-2)/2}\exp(-{\frac {1}{2}}\chi ^{2})\ \ 0\leq \chi ^{2}<\infty$

y se dice que es una variante ji-cuadrada con $n$ grados de libertad. Debido a que comenzamos con variantes normales estśndar, la densidad no contiene parámetros de ubicación o de escala sino sólo el parámetro $n$. Aún más, si las $X_{i}$ se encuentran independientemente distribuídas como variantes $N(\mu _{i},\sigma _{i}^{2})$, entonces la cantidad

$\sum _{i=1}^{n}{\dfrac {(X_{i}-\mu _{i})^{2}}{\sigma _{i}^{2}}}$

tiene una distribucion ji-cuadrada con $n$ grados de libertad, ya que cada término cuadrado en la suma ha sido normalizado con su media y su varianza. En el apéndice mostraremos una pequeña rutina en lenguaje ${\textsf {R}}$ para obtener una tabla de los valores de la distribución ji-cuadrada con varios grados de libertad. Expresaremos el porcentaje 100α de la distribución ji-cuadrada con $n$ grados de libertad como:

$\chi _{\alpha ;n}^{2}$

donde $\alpha =P(\chi ^{2}>\chi _{\alpha ;n}^{2})$.

La distribución t

La variable aleatoria $t$ con $n$ grados de libertad se define como el cociente

$t={\dfrac {z}{\sqrt {\chi ^{2}/n}}}$

de una variante normal estándar $z$ y la raíz cuadrada de una variante ji cuadrada independiente dividida esta por sus $n$ grados de libertad. La variable $t$ es adimensional y su función de densidad

$f(t)={\dfrac {\Gamma \left[(n+1)/2\right]}{{\sqrt {\pi n}}\ \Gamma (n/2)}}{\dfrac {1}{(1+t^{2}/n)^{(n+1)/2}}}\quad -\infty <t<\infty$

depende solamente del parámetro de grados de libertad $n$. Veremos más adelante como generar valores de una distribución $t$ con $n$ grados de libertad usando el lenguaje ${\textsf {R}}$.

Se expresara el porcentaje 100α como $t_{\alpha ;n}$ o

$\alpha =P(t>t_{\alpha ;n})$

La distribución F

El cociente

$F={\dfrac {\chi _{1}^{2}/m}{\chi _{2}^{2}/n}}$

de variantes ji-cuadradas independientes divididas por sus respectivos grados de libertad, tiene la distribución $F$ o cociente de varianzas con función de densidad

$f(F)={\dfrac {\Gamma \left[(m+n)/2)\right]}{\Gamma (m/2)\Gamma (n/2)}}\left({\dfrac {m}{n}}\right)^{m/2}F^{(m-2)/2}\left(1+{\frac {m}{n}}F\right)^{-(m+n)/2}\ 0\leq F<\infty$

Se denotará el porcentaje superior 100α de la distribución $F$ con $m,n$ grados de libertad por medio de $F_{\alpha ;m,n}$:

$\alpha =P(F>F_{\alpha ;m,n})$

De la definición de la variante $F$ se deduce que los puntos porcentuales inferiores se pueden obtener de los recíprocos de los valores superiores con los grados de libertad invertidos:

$F_{1-\alpha ;m,n}={\dfrac {1}{F_{\alpha ;m,n}}}$

Una Nota Breve Sobre la Funcion $\Gamma (n)$: La función $\Gamma (n)$ se define como una extensión de la función factorial a los números complejos y reales. La relación con la función factorial es:
$\Gamma (n)=(n-1)!$
La función gamma se puede definir como una integral definida en $R[z]>0$ (La forma integral de Euler):
$\Gamma (z)\equiv \int _{0}^{\infty }t^{z-1}e^{-t}\;dt$
o también
$\Gamma (z)=\int _{0}^{\infty }e^{-t^{2}}t^{2z-1}\;dt.$

Muestras Aleatorias y Estimación

Hasta este punto sólo se ha hablado de las variables aleatorias en función de su origen. Este origen son poblaciones abstractas las cuales se encuentran especificadas por su distribución o función de densidad. De vez en cuando estas funciones se les encuentra subyaciendo en algún fenómeno aleatorio, y entonces se puede crear una descripción del fenómeno a partir de su modelo matemático. Por lo general, lo que sucede realmente es que no se conocen ni la forma matemática de la distribuición ni sus parámetros. Entonces es necesario ir más alla de la teoría de la probabilidad y llegar al dominio de la inferencia estadística para obtener estimaciones de $F(x)$ o de sus parámetros a partir de muestras finitas de valores de la variable aleatoria. En esta sección se considerará un enfoque heurístico el cual lleva a obtener estimaciones con algunas propiedades requeridas.

Empecemos suponiendo que los valores de la variable aleatoria continua $X$ pueden ser observados y registrados. Aún cuando este requisito pudiera parecer obvio, no siempre es posible observar directamente la variable de interés. Pueden interferir varios factores tales como limitaciones del equipo de medición, o el fenómeno bajo observación no expresa la variable en forma directamente medible. De cualquier modo, supondremos, para simplificar las cosas, que nuestros datos son observaciones de una población continua. Supondremos también que la expresión matemática de la función de densidad $f(x;\theta _{1},\ldots ,\theta _{k})$ de $X$ se sabe a partir de algunas suposiciones anteriores aunque no se sepan los valores de los parámetros $\theta _{i}$.

Definamos ahora el espacio de parámetros de una densidad. Supongamos que la densidad depende de un sólo parámetro real $\theta$, como en el caso de una exponencial descendente. Entonces el espacio de parámetros de $\theta$ es la porción del eje de los reales que contiene todos los valores permitidos de $\theta$. En el caso de la exponencial descendente, el espacio de parámetros sería la parte positiva del eje de los reales, ya que los parámetros negativos anularían la propiedad de densidad, y un parámetro igual a cero convierte a la función en una constante. En el caso de una distribución normal con varianza conocida el espacio sería toda la recta de los reales. Extendiendo este razonamiento, el espacio de parámetros para una densidad de $k$ parámetros sería alguna región del espacio euclideano en $k$ dimensiones. Por ejemplo, el de una distribución normal sería la mitad superior del plano cartesiano con ejes $(\mu ,\sigma ^{2})$ (Fig.4).

Fig. 4 - Ejes del Espacio de Parámetros

Finalmente, se deben definir las unidades experimentales o unidades muestrales sobre las cuales se observaran los valores de $X$. Las unidades disponibles deben contener una coleccion homogénea con respecto a todas las características que puedan afectar los valores de la variable. Por ejemplo, si la variable aleatoria de interés es el nivel de triglicéridos en mujeres mexicanas adolescentes, entonces las unidades muestrales no deben incluír hombres, mujeres adultas, o niños de cualquier sexo. Las unidades muestrales deben ser independientes unas de las otras y no deben poseer cualidades comunes que puedan crear valores dependientes de $X$. Debido a que existe una abundante literatura sobre el tema de mediciones y muestreos no se abundará mas sobre el tema.

Supóngase ahora que se han seleccionado $N$ unidades aleatoriamente. Sus valores observados de $X$ se denotaran como $x_{1},\ldots ,x_{N}$. El primer problema es estimar los parametros $\theta _{i}$ a traves de una función apropiada de las observaciones. Tales estimaciones serán denotadas por un parámetro con capucha: ${\hat {\theta _{i}}}$. Estas cantidades serán llamadas estimaciones puntuales, ya que son valores asignados de $\theta$, a diferencia de los intervalos de estimación que se verán en la siguiente seccio}ón.

El problema de la estimación de parámetros no es trivial y lleva a una serie de consideraciones y cuestionamientos que hay que tomar en cuenta para la solucion del problema de estimación. Estas consideraciones se examinan a continuación:

• Estimadores no sesgados. El valor esperado del estimador debe ser igual al parametro. $E({\hat {\theta }})=\theta$.

• Consistencia. A medida que aumenta el tamaño de la muestra, ${\hat {\theta }}$ debe converger a $\theta$ con mayor probabilidad.

• Varianza minima. Muchas veces el estimador se escoge como aquel que posea la menor varianza de todos los estimadores no sesgados. Si no existe un estimador de varianza mínima para todos los tamaños de la muestra, es posible seleccionar entre diferentes estimadores a través del cociente de sus varianzas asintóticas. El estimador que tenga la mínima varianza relativa se le conoce como estimador eficiente.

• Suficiencia. Se dice que un estimador ${\hat {\theta }}$ es suficiente si contiene toda la información en las observaciones para estimar $\theta$. Esto es, el conocimiento directo de los valores $x_{1},\ldots ,x_{N}$ no proporcionará más información que la que contiene ${\hat {\theta }}$. En términos matemáticos, la distribución conjunta de las $x_{i}$ para un valor fijo de ${\hat {\theta }}$ no depende de $\theta$.

Algunas Observaciones Sobre el Sesgo:

El sesgo es la pesadilla de un estadístico. Es fácil de crear y muy díficil de domar. En ocasiones es hasta imposible de controlar. 
La definición estadística de sesgo es la sobreestima o subestima sistemática del verdadero valor. En términos un poco más claros, quiere decir que los resultados siempre estarán 
"movidos" hacia un lado en una cierta cantidad. Por ejemplo, la balanza donde me peso en la oficina del médico, siempre reporta un peso que es dos kilogramos más de lo correcto 
(en serio). Este desplazamiento consistente de dos kilogramos adicionales representa una sobreestima del peso real.
Lo más importante que hay que recordar con respecto al sesgo, es que hay que prevenirlo, o al menos minimizarlo. El sesgo es como la maleza en un jardín. Una vez que se propaga, 
es muy difícil eliminarla. Por eso es mejor eliminarla desde un principio.

Ahora bien, en cada análisis los criterios mencionados arriba deben ser evaluados con respecto al costo de obtener las muestras, la velocidad de procesamiento de los datos, y las consecuencias que puedan resultar de sesgos pequeños o varianzas grandes sobre el resultado de los experimentos. Como se verá más adelante, algunas veces es posible reducir el sesgo a traves de un escalamiento del estimador o una redefinición. En muestras pequeñas la pérdida de eficiencia en algunos estimadores puede ser compensada por su facilidad de cómputo.

Vamos a examinar ahora un método para generar estimaciones conocido como el "método de máxima posibilidad". Este procedimiento es heurístico pero se puede demostrar que los estimadores obtenidos poseen las propiedades deseadas. Empezamos definiendo la función de posibilidad de la muestra aleatoria de observaciones como la densidad conjunta de las variantes de las unidades muestrales evaluadas en $x_{i},\ldots ,x_{N}$:

$L(\theta _{1},\ldots ,\theta _{k})=\prod _{i=1}^{N}f(x_{i};\theta _{1},\ldots ,\theta _{k})$

La función de posibilidad es una medida relativa de la posibilidad de la muestra en particular $x_{i},\ldots ,x_{N}$. El método de máxima posibilidad demanda que los estimadores de $\theta _{i}$ sean seleccionados de tal manera que maximizen la función para una muestra en particular. Si la función de posibilidad tiene un máximo relativo, esto se puede lograr derivando la función y resolviendo las ecuaciones correspondientes. El valor cero de la derivada es sólo una condición necesaria para un máximo relativo y las condiciones suficientes dadas por las derivadas parciales de segundo orden deben verificarse también. Debido a que en muchas de las posibilidades que se encuentran en la inferencia estadística contienen términos exponenciales, se usara el logaritmo natural de la función de posibilidad:

$l(\theta _{i},\ldots ,\theta _{k})=\ln L(\theta _{1},\ldots ,\theta _{k})$

Entonces, si la función de posibilidad tiene un máximo relativo, las estimaciones de $\theta _{i}$ se pueden encontrar resolviendo el sistema de $k$ ecuaciones simultáneas

${\dfrac {\partial l(\theta _{1},\ldots ,\theta _{k})}{\partial \theta _{j}}}=0;\ \ j=1,\ldots ,k$

para encontrar el valor de los estimadores ${\hat {\theta _{1}}},\ldots ,{\hat {\theta _{k}}}$ Si las ecuaciones poseen raíces múltiples, será necesario escoger aquella que nos conduzca la máxima posibilidad.

Veamos dos ejemplos que describen los pasos a tomar para encontrar los estimadores de máxima posibilidad.

Ejemplo

Vamos a determinar los estimadores de máxima posibilidad de $a$ y $b$ para la densidad rectangular. La función de posibilidad de la muestra $x_{1},\ldots ,x_{N}$, es

$L(a,b)={\dfrac {1}{(b-a)^{N}}}$

Es evidente que $a$ no puede ser menor que la observación más pequeña, y que $b$ no puede ser mayor que la observación más grande. Si se denotan las observaciones ordenadas como $x_{(1)}\leq \cdots \leq x_{(N)}$,

$a\leq x_{(1)}\leq \cdots \leq x_{(N)}\leq b$

La función de posibilidad adquiere su valor máximo cuando $(b-a)$ adquiere su valor m+inimo lo que es consistente con el segundo grupo de desigualdades, y las estimaciones que minimizan el intervalo son

${\hat {a}}=x_{(1)}\;\;{\hat {b}}=x_{(N)}$

Ejemplo: La función de posibilidad de la muestra $x_{1},\ldots ,x_{N}$ de $N$ variables aleatorias normales independientes es

$L(\mu ,\sigma ^{2})={\dfrac {1}{(2\pi )^{N/2}(\sigma ^{2})^{N/2}}}\exp \left[-{\dfrac {1}{2\sigma ^{2}}}\sum _{i=1}^{N}(x_{i}-\mu )^{2}\right]$

tomando el logaritmo

$l(\mu ,\sigma ^{2})=-{\dfrac {N}{2}}\ln(2\pi )-{\dfrac {N}{2}}\ln \sigma ^{2}-{\dfrac {1}{2\sigma ^{2}}}\sum _{i=1}^{N}(x_{i}-\mu )^{2}$

Las derivadas parciales con respecto a $\mu$ y $\sigma ^{2}$ son

${\dfrac {\partial l(\mu ,\sigma ^{2})}{\partial \mu }}={\dfrac {1}{\sigma ^{2}}}\sum _{i=1}^{N}(x_{i}-\mu )$

${\dfrac {\partial l(\mu ,\sigma ^{2})}{\partial \sigma ^{2}}}=-{\dfrac {N}{2\sigma ^{2}}}+{\dfrac {1}{2\sigma ^{4}}}\sum _{i=1}^{N}(x_{i}-\mu )^{2}$

Igualando a cero y cancelando factores, queda el siguiente sistema de ecuaciones simultáneas

$\sum _{i=1}^{N}x_{i}-N\mu =0$

$\sum _{i=1}^{N}(x_{i}-\mu )^{2}-N\sigma ^{2}=0$

Despejar el estimador de $\mu$ primero, y luego usar ese valor para obtener la solución del estimador de $\sigma ^{2}$:

${\hat {\mu }}={\dfrac {1}{N}}\sum _{i=1}^{N}x_{i}={\bar {x}}$

${\hat {\sigma }}^{2}={\dfrac {1}{N}}\sum _{i=1}^{N}(x_{i}-{\bar {x}})^{2}$

$={\dfrac {1}{N}}\left[{\dfrac {1}{N}}\sum _{i=1}^{N}x_{i}^{2}-{\dfrac {1}{N}}(\sum x_{i})^{2}\right]$

Estas estimaciones son intuitivamente plausibles, ya que la media de la población se estima a través de la media de la muestra ${\bar {x}}$, y la estimación de la varianza es el promedio de la desviación cuadrada de la media.

Determinemos ahora si ${\hat {\mu }}$ y ${\hat {\sigma }}^{2}$ son estimadores insesgados. Reemplacemos las observaciones por las variables aleatorias $X$ y tomemos las esperanzas:

$E({\hat {\mu }})$

$={\dfrac {1}{N}}\sum _{i=1}^{N}E(X_{i})$

$={\dfrac {1}{N}}\sum _{i=1}^{N}\mu =\mu$

Esto significa que la media muestral es un estimador insesgado de $\mu$. El calcular la esperanza de ${\hat {\sigma }}^{2}$ involucra un poco más de cómputo:

$E({\hat {\sigma }}^{2})$

$={\dfrac {1}{N}}E\left[{\dfrac {1}{N}}\sum _{i=1}^{N}X_{i}^{2}-{\dfrac {1}{N}}\left(\sum X_{i}\right)^{2}\right]$

$={\dfrac {1}{N}}\left[NE(X_{i}^{2})-{\dfrac {1}{N}}E\left(\sum _{i=1}^{N}\sum _{j=1}^{N}X_{i}X_{j}\right)\right]$

$={\dfrac {1}{N}}\left[NE(X_{i}^{2})-E(X_{i}^{2})-{\dfrac {1}{N}}\sum \sum E(X_{i}X_{j})\right]$

$={\dfrac {1}{N}}\left\lbrace (N-1)E(X_{i}^{2})-{\dfrac {1}{N}}N(N-1)[E(X_{i})]^{2}\right\rbrace$

$={\dfrac {N-1}{N}}\sigma ^{2}$

Entonces ${\hat {\sigma }}^{2}$ no es un estimador insesgado de la varianza. Sin embargo, si reemplazamos el divisor $N$ en la fórmula original por $N-1$, se eliminará el sesgo y se obtendrá la conocida formula de la varianza de la muestra

$s^{2}={\dfrac {1}{N-1}}\left[\sum _{i=1}^{N}x_{i}^{2}-{\dfrac {1}{N}}\left(\sum x_{i}\right)^{2}\right].$

Por lo tanto, los estimadores de máxima posibilidad no son necesariamente insesgados, aunque a veces es posible corregir el sesgo multiplicándolos por un factor apropiado. Tales estimadores por lo general no tienen varianza mínima. Sin embargo, se puede demostrar que todos los estimadores de posibilidad máxima son consistentes y que sus distribuciones asintóticas son normales con parametros conocidos.

Pruebas de Hipótesis para los Parámetros de Poblaciones Normales

Pruebas Estadísticas de Hipótesis.

La inferencia estadística puede dividirse en dos áreas generales: 1) La estimación de funciones de distribución, los parámetros de tales funciones cuando se conoce su formulación matemática, o los parámetros de modelos creados alrededor de variables aleatorias; 2) El problema de probar la validez de hipótesis sobre las funciones de distribución y sus parámetros o los componentes de modelos matemáticos. En las sección anterior se estudió brevemente un aspecto de la estimación usando el principio de posibilidad máxima. Ahora se resumirán algunos principios esenciales de las pruebas de hipótesis, y posteriormente veremos como estas pruebas se pueden invertir para proporcionar intervalos de estimación de parámetros.

Quizás un ejemplo sencillo ayudará a comprender el problema de las pruebas de hipótesis en términos de las regiones del espacio de parámetros. Se sabe por experiencias pasadas que el promedio de calificaciones de los estudiantes de primer semestre en cierta universidad canadiense tienden a hallarse normalmente distribuído con media igual a 2.43 y varianza 0.04 (en algunos países las calificaciones van de 0 a 4 con incrementos decimales). Sin embargo, en este ultimo año escolar, los requisitos de entrada se hicieron más estrictos, y se presentó la hipótesis de que la calificiacion media de la población estudiantil entrante será más alta y la varianza más pequeña. Estas afirmaciones sobre los parámetros poblacionales se pueden resumir de la siguiente manera:

Hipótesis original:	$\mu =2.43$	$\sigma ^{2}=0.04$
Hipotesis alternativa:	$\mu >2.43$	$\sigma ^{2}<0.04$

'

Le llamaremos a la descripción original de la población como la 'hipótesis nula; ésta se denota normalmente como $H_{0}$. La hipótesis alternativa será denotada como $H_{1}$. Se supondrá también que la población se encuentra normalmente distribuída, pero por lo general se asume por omisión. $H_{0}$ se refiere al punto (2.43, 0.04) en el espacio de parámetros y se le conoce como una hipótesis sencilla. La hipótesis alternativa se refiere a la región sombreada de la figura 5; y ya que ese conjunto contiene mas de un punto, $H_{1}$ se conoce como una hipótesis compuesta. Un conjunto importante de hipótesis compuestas se forma por aquellos enunciados donde los valores de uno o más parámetros son completamente desconocidos. Por ejemplo, si la variable aleatoria $X$ se encuentra normalmente distribuída con varianza desconocida, las hipótesis

$H_{0}:\;\mu =\mu _{0}$

y

$H_{1}:\;\mu =\mu _{1}>\mu _{0}$

que sólo abarcan la media, se interpretarían como líneas verticales en el espacio de la figura 4 cortando al eje de las $\mu$ en $\mu _{0}$ y en $\mu _{1}$. Habría entonces que desarrollar pruebas estadísticas para tal hipótesis de tal forma que no se afectaran por el valor desconocido de $\sigma ^{2}$.

Fig. 5

A través del desarrollo de la Estadística, el objetivo de la teoría de inferencias estadísticas ha sido el desarrollar pruebas para la validez de estas hipótesis basadas en observaciones muestrales. Para el caso mas general, sea $X$ una variable aleatoria continua con los números reales como su dominio, y su función de densidad depende de los parámetros $\theta _{1},\ldots ,\theta _{k}$. El espacio muestral $W$ de todos los resultados posibles de $N$ observaciones en $X$ ser el espacio N-dimensional euclideano. Una muestra particular de observaciones se expresará en forma de vector como $\left[x_{1},\ldots ,x_{N}\right]$ y denotará un punto en el espacio muestral. Sean $\omega _{1}$ y $\omega _{2}$ dos regiones disjuntas en el espacio de parámetros $\Omega$, y se proponen las siguientes hipótesis sobre los parámetros de la distribución de $X$.

$H_{0}:[\theta _{1},\ldots ,\theta _{k}]\quad {\text{se encuentra contenida en }}\omega _{1}$
$H_{1}:[\theta _{1},\ldots ,\theta _{k}]\quad {\text{se encuentra contenida en }}\omega _{2}$

Basándonos en las observaciones muestrales deseamos decidir en alguna forma "óptima" cuál de las dos hipótesis es la más viable. Debido a que ciertos conjuntos de observaciones nos llevarían a favorecer $H_{0}$ sobre $H_{1}$ y otros conjuntos apoyarían la alternativa, la regla de decisión debe tener la siguiente forma:

${\text{Aceptar }}H_{0}{\text{ si }}\left[x_{1},\ldots ,x_{N}\right]{\text{cae en }}W-w$
${\text{Aceptar }}H_{1}{\text{ si }}\left[x_{1},\ldots ,x_{N}\right]{\text{cae en }}w$

donde $w$ es una parte en particular del espacio muestral llamada región crítica o región de rechazo de $H_{0}$. Si el verdadero estado de la naturaleza se encuentra descrito por $H_{0}$ o $H_{1}$, el tomador de decisiones puede incurrir en dos tipos de error al aplicar la regla de decisión. Un error del primer tipo, o error Tipo I, consiste en declarar $H_{1}$ como verdadera cuando en realidad $H_{0}$ es la correcta. Un error del segundo tipo, o error Tipo II se comete cuando se toma $H_{0}$ como verdadera cuando es $H_{1}$ la que describe el verdadero estado de la naturaleza. Resumimos los resultados de las acciones anteriores en la siguiente tabla:

	Estado de la naturaleza
Decisión	$H_{0}$ verdadera	$H_{1}$ verdadera
Se acepta $H_{0}$	Correcto	Error Tipo II
Se acepta $H_{1}$	Error Tipo I	Correcto

Las probabilidades de los errores de Tipo I y de Tipo II proporcionan medidas de la eficacia de la regla de decisión. La probabilidad $\alpha$ de un error de Tipo I es la probabilidad de que las observaciones muestrales se encuentren en la región crítica $w$ cuando $H_{0}$ es verdadera, o lo mismo que

$\alpha =P(\left[x_{1},\ldots ,x_{N}\right]\in w|H_{0}{\text{ es verdadera}})$

donde el símbolo $\in$ denota la inclusión en el conjunto $w$ y la barra vertical indica que el enunciado de probabilidad esta condicionado a la certeza de la hipótesis nula. El símbolo $\alpha$ se conoce como el tamaño de la prueba. La probabilidad de un error del Tipo II es

$\beta =P(\left[x_{1},\ldots ,x_{N}\right]\in W-w|H_{1}{\text{ es verdadera}})$

El complemento $1-\beta$ de la probabilidad del segundo error se conoce como la potencia de la prueba o la regla de decisión. Si la potencia se calcula para un continuo de valores de los parámetros, las probabilidades resultantes constituyen la función de potencia.

Si el tamaño de la muestra es fijo, los cambios en la forma de la región crítica que reduzcan $\alpha$ también aumentarán $\beta$. Al mismo tiempo, un aumento de $\alpha$ resultará en una reducción de $\beta$. El enfoque clásico de las pruebas de hipótesis requiere una prueba con $\alpha$ fijo cuya región de rechazo se elige de tal manera que minimize $\beta$, o que maximize la potencia. Si ambas hipótesis son sencillas, esto es, tienen la forma

	$\theta _{1}=\theta _{10}$			$\theta _{1}=\theta _{11}$
$H_{0}:$	$\cdots$		$H_{1}:$	$\cdots$
	$\theta _{k}=\theta _{k0}$			$\theta _{k}=\theta _{k1}$

entonces el lema de Neyman-Pearson afirma que la prueba más poderosa con tamaño $\alpha$ tendrá una región crítica definida por la siguiente regla de decisión:

${\text{Aceptar }}H_{0}{\text{ si }}\ \lambda ={\dfrac {f(x_{1};\theta _{10},\ldots ,\theta _{k0})\cdots f(x_{N};\theta _{10},\ldots ,\theta _{k0})}{f(x_{1};\theta _{11},\ldots ,\theta _{k1})\cdots f(x_{N};\theta _{11},\ldots ,\theta _{k1})}}>c$

y

${\text{Aceptar }}H_{1}{\text{ si }}\ \lambda <c$

donde $c$ es una constante tal que

$P(\lambda <c)=\alpha$

El lema define la región crítica como aquel conjunto de puntos en el espacio de muestras para el cual el cociente de posibilidad $\lambda$ es menor que $c$.

Nótese que el lema de Neyman-Pearson requiere que ambas hipótesis sean simples, esto es, que las regiones $\omega _{1}$ y $\omega _{2}$ del espacio de parámetros sean puntos. Las pruebas de hipótesis compuestas que involucren varios parámetros pueden ser construídas usando el poderoso criterio generalizado del cociente de posibilidad:

$\lambda ={\dfrac {L({\hat {\omega }})}{L({\hat {\Omega }})}}$

donde

$L({\hat {\Omega }})=f(x_{1};{\hat {\theta }}_{10},\ldots ,{\hat {\theta }}_{k0})\cdots f(x_{N};{\hat {\theta }}_{10},\ldots ,{\hat {\theta }}_{k0})$

es la función de posibilidad maximizada bajo la suposición de que

$H_{0}:\left[\theta _{1},\ldots ,\theta _{k}\right]\in \omega$

es cierta, y que $L({\hat {\Omega }})$ es la posibilidad maximizada con $\theta _{1},\ldots ,\theta _{k}$ la cual se permite tomar valores a través del todo el espacio de parámetros $\Omega$. Aceptamos $H_{0}$ si $\lambda >c$ y si no, se acepta la alternativa

$H_{1}:\left[\theta _{1},\ldots ,\theta _{k}\right]\in \Omega -\omega$

La constante $c$ se escoge de tal forma que $P(\lambda <c|H_{0}{\text{ es verdadera}})\leq \alpha$. Tal como en el caso de las hipótesis compuestas, cuando el verdadero tamaño de la prueba es realmente menor o igual a $\alpha$, entonces decimos que la prueba es de nivel $\alpha$. Se puede demostrar que si $H_{0}$ es verdadera, la variante

$\chi ^{2}=-2\ln \lambda$

tiende a encontrarse distribuída en forma ji-cuadrada a medida que aumenta el tamaño de la muestra con un número de grados de libertad igual a la diferencia de dimensiones entre el espacio de parámetros $\Omega$ y el subespacio de la hipótesis nula $\omega$, o lo que es mismo, el número de parámetros determinados por $H_{0}$.

Pruebas de la Media de una Variante Normal con Varianza Conocida

Sea $x_{1},\ldots ,x_{N}$ una muestra de observaciones independientes de la variable aleatoria con distribución $N(\mu ,\sigma ^{2})$. La varianza $\sigma ^{2}$ es conocida aunque $\mu$ no. Basándonos sobre las observaciones muestrales deseamos probar la hipótesis

$H_{0}:\mu =\mu _{0}$

de que la media poblacional tiene un cierto valor determinado $\mu _{0}$ contra la hipótesis alternativa:

$H_{1}:\mu =\mu _{1}>\mu _{0}$

que la media es de un valor mayor $\mu _{1}$. Usando el lema de Neyman-Pearson la prueba mas poderosa de tamaño $\alpha$ se basa en la estadística de prueba:

$z={\dfrac {({\bar {x}}-\mu _{0}){\sqrt {N}}}{\sigma }}$

cuya region crítica es

$z>z_{\alpha }$

donde $z_{\alpha }$ es el punto que divide el $100\alpha \%$ superior del resto de la distribucion normal unitaria $N(0,1)$. En términos de la media de las observaciones originales debemos rechazar $H_{0}$ en favor de $H_{1}$ si

${\bar {x}}>\mu _{0}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha }$

y si no, entonces aceptar la hipótesis nula. Si la hipotesis alternativa hubiera sido

$H_{1}':\mu =\mu _{1}<\mu _{0}$

entonces la misma estadística de prueba hubiera sido empleada, aunque la región crítica hubiera sido

$z<-z_{\alpha }$

o, equivalentemente

${\bar {x}}<\mu _{0}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha }$

Las pruebas anteriores y las hipótesis alternativas $H_{1},H'_{1}$ se conocen como hipótesis unilaterales, debido a que se indica claramente la dirección del cambio de $\mu _{0}$ a $\mu _{1}$. En la medida en que es posible hacer tales predicciones tomando en cuenta el contexto o por previas investigaciones, la potencia de las pruebas será considerablemente mayor que la de las pruebas bilaterales las cuales tienen la hipótesis alternativa:

$H_{1}^{''}:\mu =\mu _{1}\neq \mu _{0}$

lo cual permite mayores o menores valores alternativos de la media. La estadística de prueba sigue siendo , pero la región de rechazo para una prueba de tamaño $\alpha$ es de

$|z|>z_{\alpha /2}$

donde $z_{\alpha /2}$ es el punto que divide el $50\alpha \%$ superior del resto de la distribución normal unitaria. Equivalentemente, se rechaza $H_{0}$ si

${\bar {x}}>\mu _{0}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}\;\;\;\;{\text{ó}}\;\;\;\;{\bar {x}}<\mu _{0}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}$

son satisfechas.

Las curvas de potencia para las pruebas contra las tres alternativas pueden ser construídas a partir de tablas de la distribución normal, y ya que pueden ser encontradas en muchos textos de estadística, se dejará como ejercicio para la práctica del lenguaje ${\textsf {R}}$. En la práctica uno selecciona un valor apropiado de $\alpha$ y una o más probabilidades $\beta$ que sean tolerables. El tamaño de la muestra $N$ se escoge entonces para lograr que la $\beta$ mínima sea consistente con las limitaciones del presupuesto o el laboratorio.

Finalmente, se hace notar que la función de densidad normal es un miembro de la familia exponencial mencionada anteriormente en esta sección, y por lo tanto las pruebas de hipótesis compuestas tales como

$H_{0}:\mu =\mu _{0}\;\;\;\;H_{1}:\mu >\mu _{0}$

o

$H_{0}:\mu <\mu _{0}\;\;\;\;H_{1}:\mu >\mu _{0}$

son las más poderosas uniformemente.

Pruebas de Medias Cuando se Desconoce la Varianza

En muchas aplicaciones, tanto científicas, tecnológicas, financieras o de otro tipo, no es común que se conozca la varianza de la población por adelantado. Se logró un avance importante en la inferencia estadística cuando W.S. Gosset (en un artículo publicado bajo el seudónimo de "Student") obtuvo la distribución de la estadística de prueba para hipótesis compuestas sobre la media de una distribución normal con varianza desconocida. El criterio generalizado del cociente de posibilidades para probar la hipótesis contra la hipótesis alternativa basándose en $N$ observaciones independientes $x_{1},\ldots ,x_{N}$ con media ${\bar {x}}$ y varianza $s^{2}$ lleva a la estadística de prueba

$t={\dfrac {({\bar {x}}-\mu _{0}){\sqrt {N}}}{s}}$

Si $H_{0}$ es verdadera, $t$ tiene la distribución Student-Fisher con $N-1$ grados de libertad y rechazamos la hipótesis nula en una prueba con tamaño $\alpha$ si

$t>t_{\alpha ;N-1}$

donde $t_{\alpha ;N-1}$ es el punto que separa el $100\alpha \%$ superior de la distribución $t$ definida en la sección anterior. De la misma forma, si la hipótesis alternativa hubiera sido $H_{1}^{'}$, la región de rechazo estaria definida por

$t<-t_{\alpha ;N-1}$

y para la alternativa bilateral, \ref{H1biprime} la hipótesis nula será rechazada si

$|t|>t_{\alpha /2;N-1}$

La potencia de las pruebas que involucran la estadística $t$ se puede calcular de la rutina Pearson-Hartley con $\nu _{1}=1$ y $\nu _{2}=N-1$ grados de libertad y parametro de descentralidad

$\phi ={\dfrac {|\mu -\mu _{0}|}{\sigma }}\left({\dfrac {N}{2}}\right)^{1/2}$

Entonces se puede determinar, dado un $\alpha$ fijo, un tamanño de muestra que garantizará una probabilidad de potencia por encima de un mínimo especificado.

Ahora considérese que dos muestras independientes se han tomado de poblaciones normalmente distribuídas con una varianza común de $\sigma ^{2}$ pero con medias posiblemente distintas $\mu _{1}$ y $\mu _{2}$. Denótense las observaciones de las muestras como $x_{1},\ldots ,x_{N_{1}},y_{1},\ldots ,y_{N_{2}}.$ Entonces el criterio generalizado del cociente de posibilidades para la prueba de hipótesis

$H_{0}:\mu _{1}=\mu _{2}$

de medias poblacionales iguales contra la alternativa

$H_{1}:\mu _{1}>\mu _{2}$ 

lleva a la estadística de prueba

$t={\dfrac {{\bar {x}}-{\bar {y}}}{\sqrt {{\dfrac {\sum _{i=1}^{N_{1}}(x_{i}-{\bar {x}})^{2}+\sum _{i=1}^{N_{2}}(y_{i}-{\bar {y}})^{2}}{N_{1}+N_{2}-2}}\left({\dfrac {1}{N_{1}}}+{\dfrac {1}{N_{2}}}\right)}}}$

Para una prueba con tamaño $\alpha$ se rechaza $H_{0}$ en favor de $H_{1}$ si

$t>t_{\alpha ;N_{1}+N_{2}-2}$

Las regiones de rechazo para las otras hipótesis alternativas $H'_{1}:\mu _{1}<\mu _{2}$ y $H''_{1}:\mu _{1}\neq \mu _{2}$ son similares a \ref{Trej1} y a \ref{Trej2} de las pruebas con una sola muestra.

Intervalos de Confianza para Medias.

La investigadora que ha llevado a cabo un experimento costoso o complejo, rara vez se siente satisfecha al escuchar que sus observaciones han rechazado tal o cual hipótesis. Si los hallazgos muestran que el nuevo medicamento o el nuevo tratamiento tiene algún efecto significativo mas allá del efecto placebo o de las normas anteriores, la investigadora y la comunidad científica no solo preferirían saber la mejor estimación de la magnitud del efecto sino también un intervalo de valores razonables del parámetro de interés. Tales aseveraciones de los posibles valores se conocen como intervalos de confianza o, a diferencia de las estimaciones puntuales de la Sección 1.4, estimaciones de intervalo.

Supóngase que una muestra aleatoria de $N$ observaciones ha sido tomada de alguna población con densidad continua $f(x;\theta )$. El intervalo de confianza del $100(1-\alpha )\%$ para $\theta$ es aquel conjunto de valores

$t_{1}<\theta <t_{2}$

cuyos límites fueron calculados a partir de los puntos porcentuales de la función de distribución de la estimación ${\hat {\theta }}$ de $\theta$ tal que

$P(t_{1}\leq \theta \leq t_{2})=1-\alpha$

donde $\alpha$ generalmente toma el valor 0.05 o alguna probabilidad más pequeña. $1-\alpha$ se conoce como el coeficiente de confianza del intervalo. Es esencial que el enunciado de probabilidad sea interpretado como "la probabilidad de que el intervalo con extremos $t_{1},t_{2}$ contenga a $\theta$ es de $1-\alpha$" ya que en este contexto el parámetro no es una variable aleatoria. También hacemos notar que existe un número infinito de intervalos de confianza que satisfacen \ref{ConfiNt}; en muchos casos se eligirán $t_{1}$ y $t_{2}$ de tal manera que la longitud de $t_{2}-t_{1}$ sea la más corta.

Los intervalos de confianza para la media de una población normal o la diferencia de las medias de dos poblaciones normales se puede deducir de las anteriores pruebas de hipótesis. Si se han tomado $N$ observaciones independientes con media ${\bar {x}}$ de la variante $N(\mu ,\sigma ^{2})$ con varianza conocida, $H_{0}:\mu =\mu _{0}$ se rechaza en favor de

$H_{1}:\mu =\mu _{1}>\mu _{0}$

usando la prueba de tamaño $\alpha _{1}$ cuando el valor de $\mu _{0}$ está dado por

$\mu _{0}={\bar {x}}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha _{1}}$

De la misma forma, el valor más grande de $\mu _{0}$ para el cual, usando una prueba con tamaño $\alpha -\alpha _{1}$, $H_{0}$ es rechazada en favor de $H'_{1}$ cuando $H'_{1}:\mu =\mu _{1}<\mu _{0}$, entonces ese valor de $\mu _{0}$ será

$\mu _{0}={\bar {x}}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha -\alpha _{1}}$

Donde $z_{\alpha _{1}}$y $z_{\alpha -\alpha _{1}}$ son los porcentajes superiores $100\alpha _{1}$ y $100(\alpha -\alpha _{1})$ de la distribución normal unitaria. Entonces el intervalo de confianza del $100(1-\alpha )$ por ciento para $\mu$ es

${\bar {x}}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha _{1}}\leq \mu \leq {\bar {x}}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha -\alpha _{1}}$

y su longitud es de $(\sigma /{\sqrt {N}})(z_{\alpha _{1}}+z_{\alpha -\alpha _{1}})$. Pero se puede demostrar que la longitud mínima se obtiene cuando

$z_{\alpha _{1}}=z_{\alpha -\alpha _{1}}$ 

o si

$\alpha _{1}={\frac {1}{2}}\alpha$. 

Entonces el intervalo más corto con coeficiente $1-\alpha$ tiene la forma simétrica

${\bar {x}}-{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}\leq \mu \leq {\bar {x}}+{\dfrac {\sigma }{\sqrt {N}}}z_{\alpha /2}$

De la misma forma, si $\sigma ^{2}$ no es conocida, el intervalo de confianza del $100(1-\alpha )$ por ciento para $\mu$ está dado por

${\bar {x}}-{\dfrac {s}{\sqrt {N}}}t_{\alpha /2;N-1}\leq \mu \leq {\bar {x}}+{\dfrac {s}{\sqrt {N}}}t_{\alpha /2;N-1}$

donde $s$ es la desviación estándar de la muestra, y $t_{\alpha /2;N-1}$ es el $50\alpha \%$ superior de la distribución $t$ con $N-1$ grados de libertad.

Muchas veces es necesario obtener un intervalo de confianza para la diferencia de medias de dos poblaciones normales con una varianza común aunque desconocida. Las observaciones de la primera población quizá fueron obtenidas como control para aquellas en la segunda muestra que fueron obtenidas bajo un nuevo tratamiento o condición experimental. Bajo las suposiciones de normalidad y varianza común, el intervalo de confianza del $100(1-\alpha )$ por ciento para el cambio atribuible a la condición experimental es

${\bar {x}}-{\bar {y}}-s_{{\bar {x}}-{\bar {y}}}t_{\alpha /2;N_{1}+N_{2}-2}\leq \mu _{1}-\mu _{2}\leq {\bar {x}}-{\bar {y}}+s_{{\bar {x}}-{\bar {y}}}t_{\alpha /2;N_{1}+N_{2}-2}$

donde ${\bar {x}}$ y ${\bar {y}}$ son las medias de la primera y segunda muestras compuestas de $N_{1}$ y $N_{2}$ observaciones y

$s_{{\bar {x}}-{\bar {y}}}={\sqrt {{\dfrac {\sum (x_{i}-{\bar {x}})^{2}+\sum (y_{i}-{\bar {y}})^{2}}{N_{1}+N_{2}+2}}\left({\dfrac {1}{N_{1}}}+{\dfrac {1}{N_{2}}}\right)}}$

es la estimación muestral de la desviación estándar de la diferencia de medias.

Pruebas e Intervalos de Confianza para la Varianza

Las pruebas multivariables y declaraciones de confianza de los siguientes capitulos generalmente se van a aplicar a las medias y a otros parámetros de ubicación. Sin embargo, para efectos de integridad, veremos brevemente algunas hipótesis e intervalos de estimación de la varianza de una población normalmente distribuída. Si las observaciones $x_{1},\ldots ,x_{N}$ constituyen una muestra aleatoria de $N(\mu ,\sigma ^{2})$, entonces la cantidad

${\dfrac {(N-1)s^{2}}{\sigma ^{2}}}={\dfrac {\sum _{i=1}^{N}(x_{i}-{\bar {x}})^{2}}{\sigma ^{2}}}$ 

se encuentra distribuida como una variante ji-cuadrada con $N-1$ grados de libertad. El criterio generalizado del cociente de posibilidad para la prueba de la hipótesis

$H_{0}:\sigma ^{2}=\sigma _{0}^{2}$ 

contra la hipótesis alternativa

$H_{1}:\sigma ^{2}>\sigma _{0}^{2}$ 

especifica que la región de rechazo para una prueba de tamaño $\alpha$ es

${\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}>\chi _{\alpha ;N-1}^{2}$ 

donde $\chi _{\alpha ;N-1}^{2}$ es el $100\alpha$ porciento superior de la distribución ji-cuadrada con $N-1$ grados de libertad. De la misma forma, las regiones de rechazo para probar la hipótesis nula contra las alternativas $H'_{1}:\sigma ^{2}<\sigma _{0}^{2}$ y $H''_{1}:\sigma ^{2}\neq \sigma _{0}^{2}$ sería

${\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}<\chi _{1-\alpha ;N-1}^{2}$ 

y

${\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}<\chi _{1-\alpha ;N-1}^{2}\;\;o\;\;{\dfrac {(N-1)s^{2}}{\sigma _{0}^{2}}}>\chi _{\alpha _{2};N-1}^{2}$ 

respectivamente, donde en el segundo caso $\alpha _{1}+\alpha _{2}=\alpha$. Estrictamente hablando, $\alpha _{1}$ y $\alpha _{2}$ deben elegirse de tal forma que la prueba sea insesgada, esto es, su función de potencia no debe ser menor que su tamaño $\alpha$, pero para muchas aplicaciones con muestras grandes es válido partir $\alpha$ en dos partes iguales. Los intervalos de confianza para $\sigma ^{2}$ se pueden obtener directamente de la región de rechazo \ref{RejChi} como se demostrará más adelante.

Si las muestras independientes de $N_{1}$ y $N_{2}$ observaciones se han tomado aleatoriamente de las poblaciones $N(\mu _{1},\sigma _{1}^{2})$ y $N(\mu _{2},\sigma _{2}^{2})$, entonces la hipótesis

$H_{0}:\sigma ^{2}=\sigma _{0}^{2}$ 

se puede probar contra la hipótesis alternativa

$H_{1}:\sigma ^{2}>\sigma _{0}^{2}$ 

por medio de la estadística

$F={\dfrac {s_{1}^{2}}{s_{2}^{2}}}$ 

Si la hipótesis nula es verdadera, la estadística tiene una distribución $F$ con $N_{1}-1,N_{2}-1$ grados de libertad. La región de rechazo para una prueba de tamano $\alpha$ es

$F>F_{\alpha ;N_{1}-1,N_{2}-1}$ 

Por otro lado, $H_{0}$ se puede probar contra la otra hipótesis alternativa $H'_{1}:\sigma _{1}^{2}<\sigma _{2}^{2}$ por medio de la estadística

$F={\dfrac {s_{2}^{2}}{s_{1}^{2}}}$ 

cuya region crítica es

$F>F_{\alpha _{1};N_{1}-1,N_{2}-1}\ \ o'\ \ F<{\dfrac {1}{F_{\alpha _{2};N_{1}-1,N_{2}-1}}}$ 

donde la repartición $\alpha _{1}=\alpha _{2}={\frac {1}{2}}\alpha$ será suficiente.

Pruebas de Igualdad de Medias: Análisis de Varianza

Presentaremos ahora el ejemplo más sencillo del modelo lineal general que es la base del diseño experimental estadístico. La teoría de estimación y pruebas de hipótesis en modelos lineales de una sola variable ha sido discutida ampliamente en la literatura y mientras que no se presentará aquí en su totalidad, algunas partes del tema son de particular importancia para el trasfondo teórico de esta sección y servirán como una introducción al tratamiento de modelos multivariables en el capítulo correspondiente

Supóngase que se sabe que cierto compuesto bioquímico es absorbido por el cerebro aunque existe cierta evidencia que la cantidad de la sustancia por gramo de tejido cerebral parece diferir entre las cinco cepas de ratones usados en un laboratorio. Se supondrá que las cantidades relativas extraídas de los ratones sacrificados se encuentran normalmente distribuídas con la misma varianza $\sigma ^{2}$ en cada cepa. Sea $\mu _{j}$ la media poblacional para la j-ésima cepa. Es posible entonces construir una prueba de

$H_{0}:\mu _{1}=\cdots =\mu _{5}$ 

contra la hipótesis alternativa de que algunas medias son distintas usando el criterio generalizado del cociente de posibilidad. Aun mas, si $H_{0}$ es rechazada, se pueden obtener métodos para llevar a cabo pruebas simultáneas sobre las diferencias entre las medias con una probabilidad fija de error Tipo I para todas las comparaciones.

En el caso general de $k$ cepas, tratamientos, categorías de diagnóstico, o condiciones experimentales empezamos postulando que la i-ésima observación del j-ésimo tratamiento (tratamiento es el nombre que en general se le da a la característica que distingue a los $k$ grupos) se puede expresar por medio del modelo matemático

$x_{ij}=\mu +\tau _{j}+e_{ij}\ \ i=1,\ldots ,N_{j},j=1,\ldots ,k.$ 

donde

${\begin{aligned}\mu &={\text{parámetro de ubicación común a todas las observaciones}}\\\tau _{j}&={\text{efecto del j-ésimo tratamiento}}\\e_{ij}&={\text{variable aleatoria normalmente distribuída con media igual a cero y varianza}}\;\;\sigma ^{2}\end{aligned}}$

Los términos variantes $e_{ij}$ se encuentran distribuídos independientemente unos de otros. El modelo se dice que es lineal o aditivo, ya que la aplicación del j-ésimo tratamiento aumenta $\mu$ en la cantidad $\tau _{j}$, y la perturbación aleatoria también agrega o sustrae algo de los parámetros. Hacemos notar que, en una versión sencilla del ejemplo inicial, $\mu _{j}\equiv \mu +\tau _{j}$. La hipótesis de los efectos equivalentes del tratamiento se puede expresar como

$H_{0}:\tau _{1}=\cdots =\tau _{k}$ 

y se tomará como la alternativa a $H_{0}$ el modelo general \ref{GenModel} para las observaciones.

Este modelo se dice que tiene efectos fijos, ya que los $\mu _{j}$ son parámetros y cualquier inferencia a partir de las observaciones sólo se puede hacer con respecto a los $k$ tratamientos en el estudio. Si los tratamientos constituyen una muestra aleatoria de una población mayor (técnicos de laboratorio, experimentos realizados en distintos dias, etc.) entonces los efectos del tratamiento serían una variable aleatoria en sí mismos y se tendría que tomar un enfoque distinto para llevar a cabo el análisis. Esta distinción entre los modelos fijos y aleatorios de análisis de varianza fue hecha por primera vez por Eisenhart (1947) y ha sido desarrollada extensamente por otros autores en el area de diseño de experimentos. Más adelante se usaran técnicas multivariables para el análisis exacto del modelo mixto para un número fijo de tratamientos aplicados repetidamente a cada miembro de una muestra aleatoria de unidades experimentales.

Las observaciones se pueden ordenar como se muestra en la siguiente tabla

$x_{11}$	$\cdots$	$x_{1k}$
$\cdots$	$\cdots$	$\cdots$
$x_{N_{1}1}$	$\cdots$	$x_{N_{k}k}$

Denótese el total de las observaciones para el j-ësimo tratamiento como

$T_{j}=\sum _{i=1}^{N_{j}}x_{ij}$ 

y se denotará la media de ese tratamiento como ${\bar {x}}_{j}=T_{j}/N_{j}$. La suma de todas las observaciones será denotada por

$G=T_{1}+\cdots +T_{k}$ 

el número total de unidades experimentales por

$N=N_{1}+\cdots +N_{k}$ 

y la media general por ${\bar {x}}=G/N$. Expresamos entonces la suma total de los cuadrados

${\begin{aligned}S&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}(x_{ij}-{\bar {x}})^{2}\\&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}x_{ij}^{2}-{\dfrac {G^{2}}{N}}\end{aligned}}$

como la suma dos componentes independientes

$S=SST+SSE$

donde

${\begin{aligned}SST&=\sum _{j=1}^{k}N_{j}\left({\dfrac {T_{j}}{N_{j}}}-{\bar {x}}\right)^{2}\\&=\sum _{j=1}^{k}{\dfrac {T_{j}^{2}}{N_{j}}}-{\dfrac {G^{2}}{N}}\end{aligned}}$

y

${\begin{aligned}SSE&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}(x_{ij}-{\bar {x}}_{j})^{2}\\&=\sum _{j=1}^{k}\sum _{i=1}^{N_{j}}x_{ij}^{2}-\sum _{j=1}^{k}{\dfrac {T_{j}^{2}}{N_{j}}}\end{aligned}}$

Fuente	Suma de cuadrados	Grados de Libertad	Cuadrado Medio
Tratamiento	SST	$k-1$	${\dfrac {SST}{k-1}}$
Intra-tratamiento	SSE	$N-k$	${\dfrac {SSE}{N-k}}$
Total	S	$N-1$

Estos componentes se pueden resumir en la tabla de análisis de varianza mostrada en la tabla de arriba. La estadística para la prueba del cociente de posibilidad de $H_{0}$ es

$F={\dfrac {N-k}{k-1}}{\dfrac {SSE}{SST}}$ 

y, cuando $H_{0}$ es verdadera, $SST/(k-1)\sigma ^{2}$, $SSE/(N-k)\sigma ^{2}$ son variantes independientes con distribucion ji-cuadrada con $k-1$ y $N-k$ grados de libertad, respectivamente. Entonces la estadística de prueba tiene una distribución $F$ con $k-1$ y $N-k$ grados de libertad. Rechazamos $H_{0}$ con una prueba de nivel $\alpha$ si

$F>F_{\alpha ;k-1,N-k}$ 

La función de potencia para esta prueba de análisis de varianza se puede calcular a partir de la distribución F no central Pearson-Hartley. Los grados de libertad son $\nu _{1}=k-1$ y $\nu _{2}=N-k$, y el parámetro de no centralidad que mide la desviación de la hipótesis nula \ref{H0detau} de las medias poblacionales de los tratamientos es

$\phi ={\dfrac {\sqrt {\sum _{j=1}^{k}N_{j}\left(\tau _{j}-{\dfrac {\sum N_{j}\tau _{j}}{N}}\right)^{2}}}{\sigma {\sqrt {k}}}}$ 

o, como en el caso de muestras con tratamientos equivalentes que lo que generalmente se encuentra en problemas de diseño experimental

$\phi ={\dfrac {\sqrt {N_{j}\sum _{j=1}^{k}\tau _{j}^{2}}}{\sigma {\sqrt {k}}}}$ 

En esta última expresión, el término de corrección ha desaparecido de la habitual restricción $\tau _{1}+\ldots +\tau _{k}=0$ que se impone sobre los efectos del tratamiento. El parámetro $\phi$ y las gráficas de Pearson-Hartley se usan normalmente para seleccionar tamaños de muestras.

Comparaciones de Tratamientos Múltiples.

Un análisis de varianza que ha resultado en el rechazo de la hipótesis de efectos equivalentes del tratamiento no indica cuáles de esos efectos pueden ser iguales o cuales son diferentes. Este es el problema de comparaciones múltiples, o inferencias simultáneas sobre los miembros de alguna familia de hipótesis. La pruebas se construyen de tal forma que la probabilidad de error Tipo I para toda la familia será a lo mucho $\alpha$. Dos métodos serán necesarios en los siguientes capítulos, y los presentaremos a continuación en el contexto del análisis de varianza de una vía.

La primera técnica se debe a Scheffé. Defínase un contraste de parámetros llamado $\tau _{j}$ del modelo de una vía como cualquier función lineal

$\sum _{j=1}^{k}c_{j}\tau _{j}$ 

cuyos coeficientes tienen la propiedad de

$\sum _{j=1}^{k}c_{j}=0$ 

Por lo tanto $\tau _{1}-\tau _{2}$ y $3\tau _{1}+(\tau _{2}+\tau _{3}+\tau _{4})$ son contrastes mientras que $\tau _{2}-(\tau _{3}+\tau _{4})$ no lo son. En el caso particular del modelo de análisis de varianza de una vía, Scheffé ha demostrado que los intervalos simultáneos de confianza con coeficiente conjunto $1-\alpha$ para todos los contrastes $\tau _{j}$ tienen la forma

$\sum c_{j}{\bar {x}}_{j}-s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}\leq \sum c_{j}\tau _{j}\leq \sum c_{j}{\bar {x}}_{j}+s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}$ 

donde todas las sumas son sobre los $k$ tratamientos. $\sum c_{j}{\bar {x}}_{j}$ es la estimación muestral del contraste $\sum c_{j}\tau _{j}$ y

$s={\sqrt {\dfrac {SSE}{N-k}}}$ 

Nótese que $s^{2}+\sum c_{j}^{2}N_{j}$ es la estimación de la varianza del contraste estimado. Aceptamos la hipótesis nula

$H_{0}:\sum _{j=1}^{k}c_{j}\tau _{j}=0$ 

al nivel $\alpha$ si el intervalo simultáneo de confianza para ese contraste incluye el valor cero. Si, por otro lado

$\sum c_{j}{\bar {x}}_{j}>s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}$ 

o

$\sum c_{j}{\bar {x}}_{j}<-s{\sqrt {(k-1)F_{\alpha ;k-1,N-k}\sum {\dfrac {c_{j}^{2}}{N_{j}}}}}$ 

rechazamos $H_{0}$ en favor de la hipótesis alternativa unilateral $H_{1}:\sum c_{j}\tau _{j}>0$ o $H'_{1}:\sum c_{j}\tau _{j}<0$. El nivel conjunto de todas las pruebas es de $\alpha$.

La segunda técnica de comparación múltiple se conoce como el método de Bonferroni ya que se encuentra basada en una desigualdad con ese nombre. Para este procedimiento empezamos enfocándonos en una familia de $m$ declaraciones de intervalos de confianza $H_{1},\ldots ,H_{m}$ de los parámetros del modelo lineal. La probabilidad de que $H_{i}$ sea una declaración verdadera (esto es, que $H_{i}$ contenga el valor de la i-ésima función paramétrica) es igual a $P(H_{i})$ y la probabilidad de que todas las declaraciones sean correctas simultáneamente es de $P(H_{1}\cap \cdots \cap H_{m})$, donde la notación de intersección $H_{i}\cap H_{j}$ denota el evento donde ambos $H_{i}$ y $H_{j}$ son verdaderos. Esta probabilidad es el coeficiente simultáneo de confianza para la familia completa de intervalos, y nos gustaría que al menos fuera igual a algún valor especificado $1-\alpha$, donde $\alpha$ será conocida como la proporción de error para la familia de declaraciones. El cálculo de la probabilidad conjunta muchas veces es difícil para problemas prácticos y su valor depende de algunos "parámetros de molestia" que miden las intercorrelaciones de las declaraciones $H_{i}$. En lugar de calcular esta probabiliad, debemos contentarnos con la siguiente cota inferior:

$P(H_{1}\cap \cdots H_{m})\geq 1-\sum _{i=1}^{m}P({\bar {H}}_{i}).$ 

donde $P({\bar {H}}_{i})=1-P(H_{i})$ es la probabilidad de que la i-ésima declaración individual no sea verdadera. La cota es un ejemplo sencillo de las desigualdades de Bonferroni. Para un conjunto de $m$ intervalos de confianza simultáneos, por lo general se le asigna a cada declaración una proporción de error de $\alpha$, de tal manera que el coeficiente para toda la familia es de al menos $1-\alpha$.

Los intervalos de confianza de Bonferroni sobre $m$ contrastes dados

$\psi _{i}=\sum _{j=1}^{k}c_{ij}\tau _{j}\;\;\;i=1,\ldots ,m$ 

de los parámetros del modelo lineal de una vía son

$\sum _{j=1}^{k}c_{ij}{\bar {x}}_{j}-t_{\alpha /(2m);N-k}s{\sqrt {\sum _{j=1}^{k}{\dfrac {c_{ij}^{2}}{N_{j}}}}}\leq \Psi _{i}\leq \sum _{j=1}^{k}c_{ij}{\bar {x}}_{j}+t_{\alpha /(2m);N-k}s{\sqrt {\sum _{j=1}^{k}{\dfrac {c_{ij}^{2}}{N_{j}}}}}$ 

La regla de decisión para probar $H_{0}:\Psi _{i}=0$ es equivalente a aquella para el metodo de Scheffe: si (\ref{Bonf01}) no incluye el valor cero, se rechaza $H_{0}$. Para todas las comparaciones apareadas de los tratamientos, $m={\tfrac {1}{2}}k(k-1)$, mientras que para comparaciones sucesivas (suponiendo que hubo algun ordenamiento a priori), $m=k-1$. Si $m$ es pequeña los intervalos de Bonferroni pueden ser en promedio mas estrechos que los de la técnica de Scheffé, aun cuando el verdadero coeficiente de confianza de la familia es mayor que su valor nominal $1-\alpha$. Para valores de $m$ muy grandes, los intervalos de Bonferroni son imprácticamente grandes.

Ejemplo 3: En una evaluación preliminar de tres medicamentos antidepresivos, las condiciones del experimento decretaron que cada sujeto participante podía recibir solo un medicamento. Dieciocho pacientes con diagnósticos similares fueron calificados en su nivel de angustia usando una escala de cero a siete y se tomaron tres grupos de seis aleatoriamente. A cada grupo se le administró uno de los tres medicamentos, y despues de tres días los pacientes fueron calificados de nuevo usando la misma escala. Estos fueron los cambios en los niveles de angustia:

	A	B	C
	4	0	1
	5	2	0
	3	1	0
	3	2	1
	4	2	2
	2	2	2
Media	3.5	1.5	1.0

Las sumas pertinentes son

$T_{1}=21\;\;\;T_{2}=9\;\;\;T_{3}=6\;\;\;G=36\;\;\;\sum _{j=1}^{3}\sum _{i=1}^{6}x_{ij}^{2}=106$

La suma de cuadrados debida a los medicamentos es

$SST={\frac {1}{6}}(21^{2}+9^{2}+6^{2})-{\dfrac {36^{2}}{18}}=21$

y la suma total de cuadrados es $S=106-72=34$. Los valores se resumen en la tabla de análisis de varianza:

Fuente	Suma de cuadrados	Grados de libertad	Media cuadrada	F
Medicamentos	21	2	10.5	12.1
Entre Medicamentos	12	15	0.867
Total	34	17

Ya que $F_{0.01;2,15}=6.36$ concluímos que la hip[otesis de efectos equivalentes del medicamento no es válida al nivel del uno por ciento. Consultando las tablas de la distribución $F$ observamos que la $F$ excede el valor crítico para $\alpha =0.001$

Ahora usaremos el procedimiento de comparaciones múltiples de Scheffé para determinar cuáles medicamentos son los distintos. Se elige un coeficiente simultáneo de confianza de 0.99 y asi:

${\sqrt {(k-1)F_{0.01;2.15}}}=3.567$

Primero compárense los medicamentos B y C calculando el intervalo de confianza para $\tau _{2}-\tau _{3}$. Aqui $c_{1}=0,c_{2}=1,c_{3}=-1$ y la desviación estándar estimada de la población ${\bar {x_{2}}}-{\bar {x_{3}}}$ es de

${\sqrt {0.867({\frac {1}{6}}+{\frac {1}{6}})}}=0.5376$

El intervalo de confianza es $-1.42\leq \tau _{1}-\tau _{2}\leq 2.42$ y ya que esto cubre el valor cero con facilidad, podemos concluír que la hipótesis de efectos equivalentes de los medicamentos B y C no puede ser rechazada al nivel 0.01 de confianza.

De la misma forma, el intervalo simultáneo de confianza para el efecto de la diferencia entre A y B es

$0.08\leq \tau _{1}-\tau _{2}\leq 3.92$

y es posible aceptar la hipótesis alternativa

$H_{1}:\tau _{1}>\tau _{2}$

al nivel 0.01. También será interesante determinar si el medicamento A es distinto del efecto promedio de los medicamentos B y C. Aqui $c_{1}=1,\ c_{2}=-{\frac {1}{2}},\ c_{3}=-{\frac {1}{2}}$ y la desviación estándar estimada de ese contraste es

${\sqrt {0.867({\frac {1}{6}}+{\frac {1}{6}}+{\frac {1}{24}})}}=0.4653$

La estimación del contraste es de 2.25, y el intervalo de confianza es de

$0.59\leq \tau _{1}-{\frac {1}{2}}(\tau _{2}+\tau _{3})\leq 3.91$

El medicamento A parece ser distinto de los demas medicamentos restantes que acabaron siendo esencialmente los mismos.

Los intervalos de confianza de Bonferroni para los contrastes apareados con coeficiente de confianza del 0.99 son

0.13	$\leq \tau _{1}-\tau _{2}\leq$	3.87
-1.37	$\leq \tau _{2}-\tau _{3}\leq$	2.37
0.63	$\leq \tau _{1}-\tau _{3}\leq$	4.37

Introducción al lenguaje R

¿Qué es R?

El sistema ${\textsf {R}}$ para cómputos estadísticos es un ambiente para análisis de datos y gráficos. El origen de ${\textsf {R}}$ es el lenguaje ${\textsf {S}}$ el cual fue desarrollado por John Chambers y sus colegas en los laboratorios Bell por la década de los 1960s. El lenguaje ${\textsf {S}}$ fue diseñado y desarrollado como un lenguaje de programación para análisis de datos pero es, en su implementación presente, un lenguaje completo de programación con todas las características.

La distribución básica de ${\textsf {R}}$ es mantenida por un pequeño grupo de estadísticos, quienes forman el grupo central de desarrollo. Asimismo, un inmenso grupo de voluntarios le agregan al paquete una gran cantidad de funciones que aumentan enormemente su utilidad. La fuente principal de información sobre el sistema ${\textsf {R}}$ se encuentra localizada en el siguiente URL:

http://www.R-project.org

Todos los recursos de ${\textsf {R}}$ se encuentran en esta página: el sistema mismo, los paquetes accesorios, los manuales y la documentación.

La intención de este capítulo es presentar informalmente los conceptos básicos y las técnicas de manipulación de datos del lenguaje ${\textsf {R}}$ para los principiantes. No entraremos en muchos detalles técnicos. Más bien, nos concentraremos en ejemplos prácticos para que el lector pueda empezar a resolver problemas de inmediato.

Instalación del lenguaje

El sistema ${\textsf {R}}$ consiste de dos partes principales: el sistema base y la colección de paquetes accesorios que han sido contribuídos por los usuarios. El lenguaje se encuentra implementado en el sistema base. Las implementaciones de los procedimientos estadísticos y gráficos están separados del sistema base y se encuentran organizados en forma de "paquetes". Un paquete es una colección de funciones, ejemplos y documentación. La función de un paquete muchas veces se concentra en un método estadístico en particular. Ambos el sistema base y la paquetería se distribuyen a traves del "Comprehensive ${\textsf {R}}$ Archive Network" (CRAN), el cual se encuentra en

http://CRAN.R-project.org

El Sistema Base y Primeros Pasos

El sistema base se encuentra disponible en forma de código fuente y en forma binaria para varios sistemas Unix, Linux, Windows y Mac OS X. Para el analista de datos, será suficiente bajar la distribución binaria e instalarla localmente. Los usuarios de Windows pueden seguir el enlace:

http://CRAN.R-project.org/bin/windows/base/release.htm

bajar el archivo ejecutable, correrlo localmente en su máquina, y seguir las instrucciones proporcionadas por el instalador automático. Dependiendo del sistema operativo, ${\textsf {R}}$ puede ejecutarse ya sea escribiendo "R" en una terminal (Unix, Linux), o haciendo doble click en el logotipo del sistema (ver figura). %TODO: Poner logo del sistema R Para verificar que el sistema fue instalado correctamente, teclear el siguiente mando en una terminal de mandos (no incluír el signo de pesos ya que representa el apronte del sistema operativo):

 $ R --version

Aparecerá entonces en la pantalla un corto mensaje y el apronte del sistema R '>' que será algo parecido a lo siguiente:

R version 2.14.1 (2011-12-22)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-redhat-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>

Lo más probable es que la versión varíe un poco como también la plataforma, dependiente ésta de donde se ejecute el programa. Este mensaje nos indica que estamos corriendo la versión 2.14 del sistema R. Si el mensaje no aparece en la pantalla, hubo algún error y no se realizó la instalación correctamente. Consulte a su administrador de sistemas.

Se puede cambiar la apariencia del apronte del sistema ${\textsf {R}}$ con el siguiente mando:

 > options(prompt = "R> ")

Usaremos de ahora en adelante el apronte R> para mostrar los ejemplos de código a través de este libro. Un símbolo + al principio de la línea indica la continuación de un mando después de un regreso de carro.

Básicamente, el sistema ${\textsf {R}}$ evalúa los mandos proporcionados en el apronte y regresa los resultados de los cálculos. El fin de un mando se indica con la tecla de regreso (Return). Casi todos los textos de introducción a ${\textsf {R}}$ empiezan con un ejemplo usando al sistema como una calculadora:

 R> x <- sqrt(25) + 2

Este mando simple le pide al interpretador que evalúe la raíz cuadrada de 25 y le agregue 2. El resultado se lo asigna a un objeto de ${\textsf {R}}$ que tiene el nombre de variable x. El operador de asignación <- afianza el valor a su lado derecho con el nombre de la variable a su lado izquierdo. El valor del objeto x puede ser inspeccionado escribiendo

 R> x
   [1] 7

esto llama implícitamente al método print

 R> print(x)
   [1] 7

Paquetería

La distribución base incluye un conjunto de paquetes de alta prioridad como gráficos, estadísticas básicas, y algunos paquetes estándar tales como modelos lineales, regresiones, y análisis de varianza. Los paquetes que no se encuentren incluídos en la distribución base se pueden instalar directamente desde el apronte de ${\textsf {R}}$. Suponiendo que el usuario se encuentra conectado al Internet, un paquete se puede instalar proporcionando el nombre del paquete a la función install.packages. Por ejemplo, si se requiriera la estimación robusta de matrices de covarianza a través de estimadores tipo "sandwich", el paquete sandwich se puede bajar e instalar a través de:

 R> install.packages("sandwich")

Las funciones del paquete se encontrarán disponibles después de afianzar el paquete por medio de

 R> library("sandwich")

Se puede obtener una lista completa de los paquetes disponibles en

http://CRAN.R-project.org/web/packages/

Nótese que los paquetes bajados a sistemas en Windows estarán precompilados antes de bajarse e instalarse. Esto, a diferencia de los paquetes en Unix/Linux donde los paquetes se compilan localmente antes de instalarse.

Ayuda y Documentación

Existen tres formas de documentación para el sistema ${\textsf {R}}$: la ayuda en línea que viene con la distribución base, manuales electrónicos, y documentación en forma de libros y publicaciones.

La ayuda en línea consiste en una colección de páginas tipo manual de Unix (man) las cuales describen las funciones y los datos que se incluyen con ${\textsf {R}}$. Una página de manual se muestra cuando llamamos a la función help:

 R> help("mean")

o, de forma más corta:

 R> ?mean

Cada página electrónica del manual consiste de una descripción general, la lista de argumentos que requiere la función junto con una descripción de cada argumento, información sobre el valor que regresa la función, y, a veces, referencias, enlaces, y ejemplos ejecutables. La función help.search es útil para buscar dentro de las páginas manuales. Por ejemplo, se puede obtener una panorámica de los temas documentados en un paquete por medio de

 R> help(package="sandwich")

Este mando nos ayuda a obtener un resumen de las funciones del paquete sandwich.

Muchas veces un paquete viene con documentación adicional que describe las funciones del paquete y muestra ejemplos del uso del mismo. Tal documento se conoce como viñeta. Por ejemplo, la viñeta del paquete sandwich se puede abrir usando

 R> vignette("sandwich", package="sandwich")

Se puede encontrar una extensa colección de documentación en

http://CRAN.r-project.org/manuals.html

Otras funciones útiles son find y apropos. La función find nos dice en que paquete se encuentra alguna función:

> find("lowess")
[1] "package:stats"

mientras que apropos regresa con una matriz conteniendo los nombres de todos los objetos que se encuentran en la búsqueda:

  > apropos("lm")
 [1] ".__C__anova.glm"      ".__C__anova.glm.null" ".__C__glm"
 [4] ".__C__glm.null"       ".__C__lm"             ".__C__mlm"
 [7] ".__C__optionalMethod" "anova.glm"            "anova.glmlist"
[10] "anova.lm"             "anova.lmlist"         "anova.mlm"
[13] "anovalist.lm"         "colMeans"             "contr.helmert"
[16] "getAllMethods"        "glm"                  "glm.control"
[19] "glm.convert"          "glm.fit"              "glm.fit.null"
[22] "glm.nb"               "glmmPQL"              "hatvalues.lm"
[25] "KalmanForecast"       "KalmanLike"           "KalmanRun"
[28] "KalmanSmooth"         "kappa.lm"             "lm"

Ejemplos Preparados y Demostraciones

Si se desea ver el ejemplo preparado, solo hay que usar la función example usando como argumento el nombre de la función que se desea. En este ejemplo, usaremos la función de modelos lineales, o lm en esta instancia:

 R> example(lm)

El sistema presentará los resultados y gráficos incluídos en el ejemplo donde se ilustra la función lm.

Para acceder a las demostraciones de algunos paquetes el usuario puede ejecutar los siguientes mandos:

 R> demo(persp)
 R> demo(Hershey)
 R> demo(graphics)
 R> demo(plotmath)

Objetos en R

El paquete básico de ${\textsf {R}}$ incluye algunos conjuntos de datos que se pueden usar para efectos de prueba, calibración o meramente didácticos. Si se desea echar un vistazo a todos los grupos de datos que ofrece ${\textsf {R}}$ basta teclear

 R> data()

y el sistema desplegará una lista de todos los diferentes conjuntos de datos incluídos.

Para ayudar a explicar las técnicas de manipulación de datos, usaremos un grupo llamado HairEyeColor que consiste en una distribución del color de cabello, color de ojos y género en 592 estudiantes provenientes de los diferentes grupos de un curso de estadística en una universidad. Para mayor informacion sobre la naturaleza de los datos, se puede teclear

 R> ?HairEyeColor

y se desplegará la información requerida. En forma resumida, este paquete de datos consiste en los resultados de una muestra de 592 estudiantes de quienes se tomaron los datos antes mencionados. Estos datos serán útiles para ilustrar algunas técnicas para el análisis de tablas de contingencia, tales como la prueba ji-cuadrada, algunos modelos lineales y rutinas de graficado usando ${\textsf {R}}$.

Empezamos con dos mandos básicos: imprimir los datos, y analizar su estructura. Para ahorrar espacio, solamente mostraremos aquí los resultados parciales de ejecutar los mandos. Pediremos al lector que siga la lección en su propia computadora para facilitar el aprendizaje. Se teclea primero

 R> print(HairEyeColor)

y se despliega el siguiente resultado

  , , Sex = Male

       Eye
Hair    Brown Blue Hazel Green
  Black    32   11    10     3
  Brown    53   50    25    15
  Red      10   10     7     7
  Blond     3   30     5     8

, , Sex = Female

       Eye
Hair    Brown Blue Hazel Green
  Black    36    9     5     2
  Brown    66   34    29    14
  Red      16    7     7     7
  Blond     4   64     5     8

El sistema nos está mostrando como están distribuídos los colores de ojos y cabello según el género del estudiante. A continuación se teclea

 R> str(HairEyeColor)

y el sistema despliega

   table [1:4, 1:4, 1:2] 32 53 10 3 11 50 10 30 10 25 ...
 - attr(*, "dimnames")=List of 3
  ..$ Hair: chr [1:4] "Black" "Brown" "Red" "Blond"
  ..$ Eye : chr [1:4] "Brown" "Blue" "Hazel" "Green"
  ..$ Sex : chr [1:2] "Male" "Female"

Este mando presenta la estructura interna de un objeto R. En este caso nos está informando que la estructura del conjunto de datos "HairEyeColor" es una tabla de tres dimensiones: $4\times 4\times 2$. Si se desea visualizar este concepto sería como imaginar dos tablas de cuatro filas con cuatro columnas cada una. Cada una de estas dimensiones tiene un atributo al cual se le asigna un nombre. Como podemos ver en el resultado los nombre de los atributos son "Hair", "Eye", y "Sex". Los atributos son los nombres de tres distintos arreglos de caracteres. Por ejemplo, el atributo "Hair" es el nombre del arreglo de cuatro elementos que contiene los nombres de los distintos colores de cabello. Al imprimir la estructura del objeto "HairEyeColor" el sistema nos indica como están ordenados los atributos de la tabla.

Marcos de Datos (Data Frames)

Veamos ahora la otra estructura interna que ${\textsf {R}}$ usa para almacenar datos. Escribamos el siguiente mando:

 R> str(USArrests)

El sistema responde con:

  'data.frame':   50 obs. of  4 variables:
 $ Murder  : num  13.2 10 8.1 8.8 9 7.9 3.3 5.9 15.4 17.4 ...
 $ Assault : int  236 263 294 190 276 204 110 238 335 211 ...
 $ UrbanPop: int  58 48 80 50 91 78 77 72 80 60 ...
 $ Rape    : num  21.2 44.5 31 19.5 40.6 38.7 11.1 15.8 31.9 25.8 ...

${\textsf {R}}$ nos está mostrando la estructura del conjunto de datos "USArrests" el cual es el número de arrestos por cada 100,000 habitantes por los crímenes de robo, homicidio y violación correspondientes a cada uno de los 50 estados de la Unión Americana en 1973. Se proporciona también el porcentaje de la población que vive en áreas urbanas. Estos datos se encuentran estructurados en la forma de un data frame o "marco de datos". Los objetos tipo dataframe son los objetos más importantes dentro de ${\textsf {R}}$. Estos objetos representan a los datos en la forma de una tabla común. Cada fila se encuentra asociada con una sola observación y cada columna representa una variable. Las dimensiones de esta tabla se pueden extraer usando la función dim:

  R> dim(USArrests)
  [1] 50  4

indicando así que el dataframe asociado al conjunto de datos "USArrests" es una tabla con 50 filas y 4 variables. Se pueden obtener también las dimensiones de la tabla usando los mandos nrow y ncolumn:

  R> nrow(USArrests)
  [1] 50
  R> ncol(USArrests)
  [1] 4

Los nombres de las variables que describen las observaciones, esto es, los nombres de las columnas del dataframe se extraen usando

  R> names(USArrests)
[1] "Murder"   "Assault"  "UrbanPop" "Rape"

Los corchetes siempre indican un subconjunto de un objeto mayor, en este caso una sola variable extraída de la tabla completa. Debido a que los dataframes tienen dos dimensiones, siendo éstas las observaciones y las variables, se requiere la coma para especificar que se desea un subconjunto de la segunda dimensión, esto es, de las variables. Por ejemplo las observaciones correspondientes a los homicidios se encuentran representadas en un vector de longitud

  R> length(USArrests[,"Murder"])
  [1] 50

Un vector es la estructura más elemental para manejar datos en ${\textsf {R}}$ y consiste en un conjunto de elementos sencillos, los cuales son todos de la misma clase. Por ejemplo, un vector simple de los números del uno al cinco se puede construír con cualquiera de los tres siguientes mandos:

R> 1:3
[1] 1 2 3

R> c(1,2,3)
[1] 1 2 3

R> seq(from=1, to=3, by=1)
[1] 1 2 3

En el caso del conjunto de datos "USArrests", los datos se encuentran almacenados en cuatro vectores. Por ejemplo, donde se encuentran los homicidios, es un vector numérico:

  R> class(USArrests[,"Murder"])
  [1] "numeric"

El tamaño del vector donde se encuentran almacenados

  R> length(USArrests[,"Murder"])
  [1] 50

Y el primer elemento del vector es

  R> USArrests[,"Murder"][1]
  [1] 13.2

Las mediciones a nivel nominal o categorías se encuentran representadas por variables de tipo factor tales como la categoría a la cual pertenece un cierto tratamiento. Por ejemplo, el conjunto de datos "InsectSprays" representa el conteo de plagas en ciertas unidades agrícolas experimentales las cuales fueron tratadas con diferentes insecticidas nombrados aquí como "A", "B", "C", "D", "E" y "F". Veamos la estructura de este conjunto de datos:

  R> str(InsectSprays)
'data.frame':   72 obs. of  2 variables:
 $ count: num  10 7 20 14 14 12 10 23 17 20 ...
 $ spray: Factor w/ 6 levels "A","B","C","D",..: 1 1 1 1 1 1 1 1 1 1 ...

Los objectos de tipo factor y los de tipo character difieren en la forma en la que son almacenados internamente. Cada elemento de un vector de clase character se almacena como una variable character mientras que un objeto de tipo factor es almacenado como una variable entera indicando el nivel del factor. En otras palabras, el nombre de la categoría usada se convierte en un número entero que va del 1 al número de posibles valores de la categoría. En el caso de los insecticidas, el valor numérico de la variable spray corre del 1 al 6.

  R> nlevels(InsectSprays[,"spray"])
  [1] 6
  R> levels(InsectSprays[,"spray"])
  [1] "A" "B" "C" "D" "E" "F"

Como ejemplo de una estadística sencilla, se puede obtener un conteo de la frecuencia de los niveles de cada variable tipo factor a partir de

  R> table(InsectSprays[,"spray"])

 A  B  C  D  E  F
12 12 12 12 12 12

Importación y Exportación de Datos

En la sección anterior vimos como leer datos en línea y usar los paquetes de datos incluídos con el lenguaje. Ahora veremos formas más relevantes de como importar datos al sistema ${\textsf {R}}$. Los formatos más frecuentes que se encuentran en la práctica son las hojas de trabajo de Microsoft Excel, archivos de texto, y conexiones con bases de datos relacionales usando SQL. El importar datos directamente de una base de datos tal como Oracle o Microsoft SQL Server implica el conocimiento de SQL y por lo tanto dejaremos el tema para un capítulo posterior. Por lo pronto nos abocaremos a las hojas de trabajo y a archivos de texto. Los archivos de texto pueden venir en dos formas: campo fijo (donde cada elemento ocupa un lugar fijo de espacios) o campo variable delimitado por un carácter. Estos últimos por lo general vienen delimitados por una coma, un espacio, un tabulador, o un par de comillas. Supongamos que tenemos un archivo de texto el cual contiene información que queremos procesar. El archivo lleva el nombre de users.csv y es un archivo que contiene los datos de un grupo de usuarios. El archivo sólo contiene un identificador numérico, la edad del usuario, género, profesión, y el código postal de su domicilio. Los elementos se encuentran separados por comas y los datos de tipo texto se encuentran delimitados por comillas. Mostramos algunos renglones del archivo para obtener una idea de su estructura:

880,13,"M","student","83702"
881,39,"M","marketing","43017"
882,35,"M","engineer","40503"
883,49,"M","librarian","50266"
884,44,"M","engineer","55337"

Usaremos entonces la función read.table para leer el archivo:

R> users <- read.table("users.csv", header=TRUE, sep=",")

El argumento header=TRUE indica que la primera línea del archivo users.csv debe interpretarse como los nombres de las variables a leerse. Las columnas se encuentran separadas por una coma (sep=",") y los datos se almacenarán en el dataframe users. Se puede usar también la función read.csv para leer archivos separados por comas pero la función read.table es un poco más general y es mejor que nos acostumbremos a usarla. La función read.table trata de adivinar cuál es la clase de cada variable al leer el archivo. Veamos como fue estructurado dentro de ${\textsf {R}}$:

  R> str(users)
'data.frame':   943 obs. of  5 variables:
 $ UserID    : int  1 2 3 4 5 6 7 8 9 10 ...
 $ Age       : int  24 53 23 24 33 42 57 36 29 53 ...
 $ Gender    : Factor w/ 2 levels "F","M": 2 1 2 2 1 2 2 2 2 2 ...
 $ Occupation: Factor w/ 21 levels "administrator",..: 20 14 21 20 14 7 1 1 19 10 ...
 $ Zipcode   : Factor w/ 795 levels "00000","01002",..: 623 690 271 332 134 759 652 49 2 643 ...

Podemos observar que cada observación (cada línea del archivo) está compuesta por cinco variables: dos de clase int (enteros) y tres variables de tipo factor cada una con distinto número de niveles.

La exportación de datos desde ${\textsf {R}}$ se hace de forma similar. Se construye un archivo delimitado con comas para que pueda ser leído por Microsoft Excel a partir de un objeto tipo dataframe a través de:

  R> write.table(users, file="USERS.csv", sep=",", col.names=NA)

Si se desea guardar los datos para que solo sean procesados por ${\textsf {R}}$, cualquier objeto en ${\textsf {R}}$ se puede salvar en un archivo binario externo a través de:

  R> save(users, file="USERS.rda")

donde la extensión .rda es estándar. Si deseamos cargar los objetos de un archivo tipo rda a ${\textsf {R}}$.

 R> load("USERS.rda")

Manejo Básico de Datos

Los ejemplos mostrados en la sección anterior muestran la importancia de usar los dataframes para el manejo y almacenamiento de datos tabulares en ${\textsf {R}}$. Internamente, un dataframe es una lista de vectores con una longitud común $n$, donde $n$ es el número de filas en la tabla. Cada uno de estos vectores representa las medidas de una variable y hemos vista que se puede acceder tal variable por medio de su nombre. Por ejemplo, los tipos de profesiones en el juego de datos users.csv:

  R> occ <- users[,"Occupation"]
  R> occ[1]
  [1] technician

El vector occ es de clase Factor con 21 niveles y tiene una longitud de 943 elementos. Se puede extraer un subconjunto de los elementos del vector occ usando el operador []. Por ejemplo, para extraer los elementos 125 al 128, usamos el siguiente operador:

  R> occ[125:128]
  [1] lawyer    lawyer    none      marketing
  21 Levels: administrator artist doctor educator engineer ... writer

${\textsf {R}}$ nos recuerda que la variable que estamos leyendo es de tipo Factor y que posee 21 niveles que van desde "administrator" a "writer". Nota: El llamarles "niveles" a los posibles valores que puede tomar la variable Occupation no implica que exista una jerarquía entre ellos. Es simplemente un asunto de nomenclatura.

Una característica interesante de los índices en ${\textsf {R}}$ es el uso de índices negativos. El uso de éstos indica todos los elementos que no son parte del vector indicado por los corchetes. Por ejemplo, si queremos indicar los primeros cinco elementos del vector occ' lo podemos hacer de la siguiente forma usando índices negativos:

  R> occ[-(6:943)]
  [1] technician other      writer     technician other

Y en forma convencional con índices positivos:

  R> occ[1:5]
  [1] technician other      writer     technician other

Se puede imprimir una tabla de los primeros tres registros (filas) del juego de datos usando el hecho de que los dataframes manejan el concepto de filas y columnas. Entonces, separamos los subconjuntos correspondientes a las filas y a las columnas por medio de una coma. Por ejemplo:

  R> users[1:3, c("UserID", "Age", "Gender", "Occupation", "Zipcode")]
    UserID Age Gender Occupation Zipcode
  1      1  24      M technician   85711
  2      2  53      F      other   94043
  3      3  23      M     writer   32067

Con esta instrucción, le estamos diciendo a ${\textsf {R}}$ que imprima los tres primeros registros (filas) de los datos con las columnas indicadas. Al mismo tiempo, si quisiéramos extraer una sola variable del dataframe usaríamos:

  R> occ <- users$Occupation
  R> occ[1:5]
  [1] technician other      writer     technician other
  21 Levels: administrator artist doctor educator engineer ... writer

La instrucción mostrada arriba es equivalente a la instrucción

  R> occ <- users[,"Occupation"]

Ejemplo: Nos interesa saber cuáles son los datos correspondientes a cinco usuarios de mínima edad. Primero, ordenaremos las edades de menor a mayor, usando los índices correspondientes del vector Age para usarlos como índices en el vector users.

  R> ages <- order(users$Age)
  R> users[ages[1:3], c("UserID", "Age", "Gender", "Occupation", "Zipcode")]
    UserID Age Gender Occupation Zipcode
30      30   7      M    student   55436
471    471  10      M    student   77459
142    142  13      M      other   48118
609    609  13      F    student   55106
289    289  11      M       none   94619

La función order ordena el vector Age de menor a mayor y regresa los índices del vector en este nuevo orden. Usamos entonces estos índices para extraer los cinco elementos de users que poseen las cinco mínimas edades.

Otra forma de extraer datos es usar un vector lógico que adquiere el valor TRUE cuando selecciona el elemento deseado y FALSE cuando no. Por ejemplo, para obtener los datos de aquellos usuarios menores de 15 años, expresaríamos:

  R>  users[users$Age<15,  c("UserID", "Age", "Gender", "Occupation", "Zipcode")]
      UserID Age Gender Occupation Zipcode
30      30   7      M    student   55436
142    142  13      M      other   48118
206    206  14      F    student   53115
289    289  11      M       none   94619
471    471  10      M    student   77459
609    609  13      F    student   55106
628    628  13      M       none   94306
674    674  13      F    student   55337
813    813  14      F    student   02136
880    880  13      M    student   83702
887    887  14      F    student   27249
 R>

donde la expresion users$Age <15 indica un vector lógico de longitud 943 donde

R> table(users$Age<15)

FALSE  TRUE
  932    11

los elementos son falsos(FALSE) or verdaderos (TRUE).

Si tuviéramos datos faltantes, éstos estarían indicados por medio de un símbolo especial: NA indicando así que esa medición no se encuentra disponible. Si queremos saber cuáles elementos faltan en alguna columna usamos la siguiente función

 R> na_occ <- is.na(users$Occupation)
 R> table(na_occ)
na_occ
FALSE
  943

La función regresó con un vector lógico de 943 elementos todos ellos igual a FALSE. Lo que quiere decir que no hay elementos faltantes en la columna Occupation.

El tomar subconjuntos de los dataframes usando expresiones lógicas nos puede ahorrar mucho tecleado. La función subset toma un dataframe como su primer argumento, y una expresión lógica como el segundo. Por ejemplo, podemos seleccionar un subconjunto de todos los usuarios mayores de 40 años por medio de

 R> mayores50 <- subset(users, Age> 50)
 R> dim(mayores50)
[1] 105   5

El vector mayores50 es un vector de 105 filas con 5 columnas. Nótese que no es necesario especificar en el primer argumento la variable Age bajo consideración.

Cálculos con Datos

Algunas Estadísticas Sencillas

Dos funciones son muy útiles para extraer información general sobre los objetos en ${\textsf {R}}$: str y summary. La primera es para obtener información sobre el tipo y cantidad de los objetos que se encuentran almacenados en el dataframe; y la segunda es para obtener estadísticas de sumarización del conjunto de datos.

  R> summary(users)
     UserID           Age        Gender          Occupation     Zipcode
 Min.   :  1.0   Min.   : 7.00   F:273   student      :196   55414  :  9
 1st Qu.:236.5   1st Qu.:25.00   M:670   other        :105   55105  :  6
 Median :472.0   Median :31.00           educator     : 95   10003  :  5
 Mean   :472.0   Mean   :34.05           administrator: 79   20009  :  5
 3rd Qu.:707.5   3rd Qu.:43.00           engineer     : 67   55337  :  5
 Max.   :943.0   Max.   :73.00           programmer   : 66   27514  :  4
                                         (Other)      :335   (Other):909

Podemos ver que para la variable numérica Age, la función arrojó los valores extremos, el primer y tercer cuartil, la mediana y la media. Para el caso de las variables nominales Gender, Occupation, y Zipcode, la función calculó el número correspondiente a los niveles de la variable. Esto es, la variable Gender tiene 273 valores de F, y 670 de M. En el caso de las variables Occupation y Zipcode, la función proporciona los 6 valores más frecuentes.

Programando con R

En las secciones anteriores hemos estado haciendo todos nuestros cálculos con las funciones internas de ${\textsf {R}}$. Sin embargo, el verdadero poder de ${\textsf {R}}$ proviene de la expresividad de este lenguaje en forma de funciones que se usan para las tareas de análisis. Estudiaremos dos ejemplos pequeños a continuación.

Digamos que queremos agregar una medida robusta de variabilidad a las medidas de ubicación calculadas en la sección anterior. Aparte de la media y la mediana, queremos calcular el intervalo intercuartil. Esto es, la diferencia entre el tercer y el primer cuartil. Haciendo una búsqueda rápida del manual en línea (usando help("interquartile")) obtenemos el resultado de IQR, pero en lugar de usar directamente esta función, usaremos la función quartile para calcular cuantiles.

Una función en ${\textsf {R}}$ no es más que un objeto, y todos los objetos se crean de la misma manera. Esto es, tenemos que asignarle un objeto de tipo function a una variable. Un objeto function consiste de una lista de argumentos, valores por omisión, y un cuerpo que define los cálculos. El cuerpo empieza y termina con llaves {} y, por supuesto, se entiende que el cuerpo consiste de código en ${\textsf {R}}$

Regresando a nuestro ejemplo, le llamaremos a nuestra función iqr. La función iqr deberá operar sobre vectores numéricos, por lo que tendrá como argumento de entrada un vector x. Este vector numérico le será pasado a la función quantile para calcular los cuartiles en la muestra. La diferencia requerida entre el tercer y el primer cuartil se calculará usando diff. La definición de nuestra función se presenta a continuación:

R> iqr <- function(x) {
       q <- quantile(x, prob = c(0.25, 0.75), names = FALSE)
	return(diff(q))
}

Una prueba usando datos simulados provenientes de una distribución normal estándar muestra que nuestra función sí opera como debido. Si la comparamos con la función IQR mostramos que el resultado es el correcto

R> xdata <- rnorm(100)
R> iqr(xdata)
[1] 0.9977416
R> IQR(xdata)
[1] 0.9977416

Sin embargo, cuando el vector numérico contiene valores faltantes, la función falla como se muestra en el siguiente ejemplo:

R> xdata[1] <- NA
R> iqr(xdata)
Error in quantile.default(x, prob = c(0.25, 0.75), names = FALSE) :
  missing values and NaN's not allowed if 'na.rm' is FALSE

Para hacer la función un poco más flexible, sería conveniente agregar todos los argumentos de quantile a la lista de argumentos de iqr. Si lo hacemos por simple copia podría llevar a inconsistencias y errores como, por ejemplo, cuando la lista de argumentos de quantile cambia. En lugar de eso, usamos el argumento de puntos, el cual sirve de comodín para cualquier argumento de la función

R> iqr <- function(x, ...) {
+   q <- quantile(x, prob = c(0.25, 0.75), names=FALSE,
+                 ...)
+   return(diff(q))
+ }
R> iqr(xdata, na.rm=TRUE)
[1] 1.019106
R> IQR(xdata, na.rm=TRUE)
[1] 1.019106

El argumento na.rm=TRUE detecta si existen valores faltantes (NA) en el vector de entrada y los extrae para efectos del cálculo. El vector de entrada no es alterado al salir de la función.

Gráficas Sencillas

El grado de oblicuidad (skewness) de una distribución se puede obtener construyendo histogramas usando la función hist (Existen alternativas un poco más sofisticadas pero se verán en capítulos posteriores). Por ejemplo, el código para crear la figura 6 primero divide la región de graficando horizontalmente en dos partes igualmente espaciadas (layout) y luego grafica el histograma de las edades en la parte superior y el histograma del logaritmo de la edad en el segundo renglón. Este último parece ser un poco más

simétrico.

Histogramas de edades

Las relaciones bivariantes de dos variables continuas se grafican generalmente como una gráfica de puntos. En ${\textsf {R}}$, las relaciones de regresión se especifican a través de fórmulas modelo los cuales, en el caso de dos variables, pueden tomar una forma parecida a

  R> fm <- ingresos ~ ventas
  R> class(fm)
   [1] "formula"

donde la variable dependiente es la izquierda y la independiente la derecha (ventas). El tilde separa las partes izquierda y derecha. Esta fórmula se puede pasar como argumento a una función modelo como parte de un cálculo más complicado.

Supongamos que ahora deseamos las distribuciones de las edades pero separadas por género. Esto es, queremos las distribuciones de las edades masculinas y las femeninas por separado. El código ${\textsf {R}}$ para obtener una gráfica de cajas de las edades es:

R> tmp <- subset(users, Gender %in% c("M", "F"))
R> plot(Age ~ Gender, data = tmp, ylab="Edades", xlab="Genero", varwidth=TRUE)

Podemos observar el resultado en la figura 7

Gráfica de cajas (Boxplots) por edades

.

Organizando un Análisis

Aún cuando es posible realizar un análisis estadístico tecleando todos los mandos directamente en el apronte de ${\textsf {R}}$, es mucho más cómodo tener un archivo de texto por separado que contenga todos los pasos necesarios para llevar a cabo el análisis requerido. El archivo, que podemos llamar análisis.R, se puede correr dentro del ambiente de ${\textsf {R}}$ usando el mando source:

  R> source("analisis.R", echo = TRUE)

Si ponemos todos los mandos juntos necesarios para un análisis, podemos reproducir el análisis a voluntad con mayor número de datos o con distintos datos según la necesidad.

Algebra Matricial

Introducción

En un capítulo posterior veremos que una variable aleatoria multidimensional es simplemente un conjunto ordenado $[X_{1},\ldots ,X_{p}]$ de variantes sencillas. Cuando decimos "ordenado" queremos decir que la variante que describe la i-ésima faceta de cada unidad muestreada de la población siempre aparecerá en la i-ésima posición de la secuencia. El número de variantes $p$ en el arreglo será constante a través del problema o el análisis que se esté llevando a cabo. Por ejemplo, supóngase que los componentes de la variante son los pesos al nacer de los ratones hembra de una cierta cepa y los demás componentes son los pesos a los 10, 20 y 30 días de nacidas. Los pesos entonces son descritos por una variable aleatoria $[X_{1},X_{2},X_{3},X_{4}]$ con alguna distribución en el espacio de 4 dimensiones. Si esos pesos fueron registrados con una muestra de $N$ ratones, las observaciones se pueden resumir por medio del arreglo

${\begin{bmatrix}x_{11}&\cdots &x_{14}\\\cdots &\cdots &\cdots \\x_{N1}&\cdots &x_{N4}\end{bmatrix}}$

donde las filas corresponden a los distintos ratones. Tales arreglos rectangulares y lineales de números se conocen como matrices y vectores respectivamente, y las reglas para su manipulación se conocen como álgebra lineal. Esta álgebra es el lenguaje virtual del análisis multivariable, especialmente en la parte más común e importante la cuál se encuentra basada en la distribución normal multidimensional. Incluso, se puede afirmar que las técnicas de análisis multivariable no hubieran podido ser desarrolladas sin la facilidad, la conveniencia y la elegancia de las matrices.

En este capítulo se resumirán algunas propiedades, operaciones y teoremas del álgebra matricial necesarias para el material que se avecina. Si se desea profundizar en el tema, se proporcionarán referencias para consulta y aprendizaje.

Algunas Definiciones

Supongamos que tenemos a nuestra disposición una serie de elementos que se comportan de acuerdo a un conjunto de axiomas, por ejemplo, los números reales o complejos. Definimos una matriz

$A={\begin{bmatrix}a_{11}&\cdots &a_{1c}\\\cdots &\cdots &\cdots \\a_{r1}&\cdots &a_{rc}\end{bmatrix}}$

como un arreglo rectangular ordenado de los elementos. El término general de la matriz se expresará como $a_{ij}$ donde el primer subíndice se refiere a la i-ésima fila y el segundo subíndice se refiere a la j-ésima columna. Las dimensiones de una matriz son importantes, de modo que una matriz con $r$ filas y $c$ columnas será referida como una matriz de $r\times c$ dimensiones.

Para diferenciarlos de las matrices, les llamaremos a los números y variables unidimensionales que manejamos habitualmente, escalares. Un escalar es considerado como una matriz de $1\times 1$. En las primeras secciones de este capítulo supondremos que los elementos de las matrices son escalares reales. Posteriormente trataremos con matrices cuyos elementos son matrices de dimensiones más pequeñas.

Un vector es una matriz con una sola fila o columna. Normalmente, expresaremos el vector columna de $n$ componentes como

$x={\begin{bmatrix}x_{1}\\\vdots \\x_{n}\end{bmatrix}}$

De la misma forma, el vector fila lo expresaremos como

${\textbf {x}}^{\prime }=[x_{1},\ldots ,x_{n}]$

y consiste de una sola fila de $n$ elementos. Cada vector especifica las coordenadas de un punto en el espacio euclideano $n$-dimensional y la conexión con el concepto físico o analítico de un vector es aparente si se considera ese punto como el término de un segmento de recta que empieza en el origen de los ejes coordenados.

El apóstrofe que acompaña al símbolo del vector fila indica que ${\textbf {x}}^{\prime }$ es la transpuesta del vector columna x. En general, si A es cualquier matriz de $r\times c$ dimensiones, la transpuesta de A es la matriz ${\textbf {A}}^{\prime }$ de $c\times r$ dimensiones formada al intercambiar las filas y las columnas:

${\textbf {A}}^{\prime }={\begin{bmatrix}a_{11}&\cdots &a_{1c}\\\cdots &\cdots &\cdots \\a_{r1}&\cdots &a_{rc}\end{bmatrix}}={\begin{bmatrix}a_{11}&\cdots &a_{r1}\\\cdots &\cdots &\cdots \\a_{1c}&\cdots &a_{cr}\end{bmatrix}}$

Si una matriz es cuadrada e igual a su transpuesta, se dice que es simétrica. Entonces $a_{ij}=a_{ji}$ para todos los pares de $i$ y