Plan de Trabajo

Overview

Teaching: 0 min
Exercises: 0 min

Questions

• Key question (FIXME)

Objectives

• First learning objective. (FIXME)

Plan de Trabajo

Este laboratorio se ofrece con el propósito de capacitar al oyente para realizar inferencias básicas sobre conjuntos de datos estadísticos con la ayuda de las tecnologías de computación disponibles en la actualidad.

Para este fin, se expondrán los conceptos fundamentales sobre los siguientes tópicos.

• Introducción al lenguaje R (2 horas).
• Distribuciones de probabilidad que dependen de uno o más parámetros (2 horas).
• Métodos de estimación de parámetros, puntual y por intervalos (2 horas).
• El concepto de prueba de hipótesis estadística y ejemplos de las pruebas más comunes. En particular las pruebas de bondad de ajuste (6 horas).
• Regresión lineal simple (6 horas).
• Análisis de varianza (2 horas).

En total el laboratorio se desarrollará en alrededor de 20 horas de trabajo en grupo.

Key Points

• First key point. Brief Answer to questions. (FIXME)

Introducción al lenguaje R

Overview

Teaching: 0 min
Exercises: 0 min

Questions

• Key question (FIXME)

Objectives

• First learning objective. (FIXME)

Comandos básicos

El lenguaje R es sensible a las mayúsculas: no es lo mismo variable1 que Variable1.

Para interrumpir un cálculo extenso, oprimir Esc.

Para salir de R:

  q()

Asigna el valor “5” a la variable $ x $ :

  x<-5

Da la lista de las variables, u objetos, definidos por el usuario:

  ls()

Elimina la variable $ x $ :

    > rm("x") 

Elimina todos los objetos definidos por el usuario:

    > rm(list=ls())

Muestra el directorio en uso:

    > getwd()

Cambia el directorio en uso:

    >setwd("/media/euba/ADATA UFD/Diplomado/Programas")

Muestra los archivos en el directorio en uso:

    > dir()

Nos da los últimos 10 comandos ejecutados:

    > history(10)

Guarda la lista de comandos ejecutados en el archivo borrar.txt:

    > savehistory("borrar.txt") 

Carga el archivo borrar.txt:

    > loadhistory("borrar.txt")

R puede ejecutar un conjunto de instrucciones leyéndolas desde un archivo externo. Por ejemplo, considere el archivo Sumar.r que se muestra a continuación:

    > #suma
    > a<-1
    > b<-2
    
    > print(a+b)

Para ejecutar el archivo Sumar.r se usa el comando:

    > source("Sumar.r")

Ambiente de Cálculo

Define el vector $ y $ :

    > y <-c(1, 2, 3, 4, 5, 6, 7, 8, 9)
    > y 

Y despliega lo siguiente:

    > [1] 1 2 3 4 5 6 7 8 9

Define la matriz de $ 3x3 $ :

    > m <-matrix(y, 3, 3) 
    > m

    >     [,1]    [,2]    [,3]
      [,1]  1       2       3
      [,2]  4       5       6
      [,3]  7       8       9

Define la matriz $ u $ de $ 3x1 $

    > u<-matrix(1, 3, 1)
    > u

    >   [,1]
     [1,] 1
     [2,] 1
     [3,] 1

Calcula la media de las entradas en el vector $ y $ :

    > mean(y)

En la siguiente tabla se reportan los comandos para calcular otros estadísticos básicos.

En la siguiente tabla se reportan algunas operaciones matemáticas.

Ejemplos

    > y <- c(1:5)

Sirve para generar el vector:

    > [1] 1 2 3 4 5

    > seq(from=1, to=13, by=3)

Da como resultado:

    > [1] 1 4 7 10 13

Si $ x $ es un vector, entonces $ x[i] $ es la i-ésima entrada de $ x $ . Por ejemplo:

    > x <-c(4, 6, 2, 4, 8, 9, 1, 5)
    > x[3]
    [1] 2
    > x[1:3]

Da como resultado:

    > [1] 4 6 2

Gráficas

Especifica las dimensiones de una gráfica

    > dev.new(width=6, height=5)

Especifica el tamano de los carácteres de la grafica:

    > par(cex=1.5)

Produce la gráfica:

    >x<-1:10
    >y<-x^2
    >plot(x, y)

La gráfica anterior se puede imprimir en un archivo pdf mediante las siguientes instrucciones:

    >pdf("rplot.pdf", width=6, height=5)
    >par(cex=1.5)
    >plot(x, y)
    >dev.off()

Las funciones para graficar en R tienen varios parámetros. El comando par sirve para especificar de manera global los parámetros de las gráficas que se realizarán en la sesión. Por ejemplo:

    >par(cex=1.5) 

Sirve para definir el tamaño de los caracteres de las gráficas. xlab y ylab sirven para especificar las leyendas de los ejes horizontal y vertical.

    >plot(x, y, xlab="equis", ylab="cuadrado")

col y bg cambian el color de la gráfica y el color del fondo de la gráfica.

    >par(bg="yellow")
    >plot(x, y, col="red")

pch para seleccionar el símbolo para los puntos de la gráfica.

    >plot(x, y, pch=1:10)

lty y lwd cambian el tipo de línea y el grueso de la línea respectivamente. lines sirve para a ̃nadir a la gráfica una curva a trazo continuo.

    >plot(x, y, lwd=3)
    >lines(x, x ˆ 2, lty=5)

cex tamaño del texto y puntos de la gráfica.

    >plot(x, y, cex=2:3)

main es el título de la gráfica.

    >plot(x, y, main="hola")

ps sirve para especificar el tamaño del texto dentro de la gráfica.

    >par(ps=10, cex=1.5, cex.main=2)
    >plot(x, y, cex=2:3, main=“Cuadrado”)

fg especifica el color del marco de la gráfica.

    >plot(x, y, fg="blue")

xlim y ylim especifica el rango en el eje de la $ x $ y en el eje de las $ y $ , respectivamente.

    >plot(x, y, ylim=c(-10, 110), xlim=c(-1, 12))

text sirve para añadir texto a la gráfica.

    >text(4, 20, "(4, 16)")

mtext añade texto en los márgenes.

    > mtext("aqui", side=4)

type especifica el tipo de gráfica.

    > plot(x, y, type = "b")

La gráfica puede ser de uno de los tipos que se especifican en la siguiente tabla.

Ejemplos Adicionales

    > plot(x, y, pch=21, lwd=2, col="red", bg="blue")

    > plot(x, y, pch=c("a", "b"), col=c("red", "blue"))

    > par(cex.lab=2, cex.main=2)

Tamaño de las leyendas en los ejes y título

    > plot(x, y, pch=c("a", "b"), col=c("red", "blue"), +main="cuadrado")

Generación de Números Aleatorios

Se inicializa el generador de números aleatorios

    > set.seed(777)

Si el generador de números aleatorios se inicia con un mismo valor cada vez que se corra el comando, entonces la sucesión de números aleatorios permanecerá sin cambiar.

    > runif(5)

Da como resultado:

    > [1] 0.6878574 0.4921926 0.3451156 0.9950499 0.6952672

R maneja distintos tipos de variables aleatorias. En la siguiente tabla se reportan los nombres de algunas de las variables aleatorias de uso común.

Los nombres de las variables aleatorias se usan en conjunto con una raíz que indica la función a ejecutar.

    > rnorm(5, 1, 0.01)

Da como resultado:

    [1] 0.9979378 0.9962103 0.9969574 1.0005416 0.9811907

    > rbinom(5, size=10, prob=0.5)

Da como resultado:

    [1] 5 5 9 6 7

La siguiente instrucción sirve para calcular la probabilidad de que una variable aleatoria $ X \sim Binomial(10, 5) $ sea igual a 7.

    > dbinom(7, size=10, prob=0.5)

    [1] 0.1171875

Nota

R reconoce a dbinom como a una función de densidad, aún cuando en sentido estricto se trata de una función de distribución de masa.

Con la siguiente instrucción se calcula la probabilidad de que X \sim Binomial(10, 5) sea menor o igual a 7.

    > pbinom(7, size=10, prob=0.5)

Da como resultado:

    [1] 0.9453125

Para calcular los cuantiles:

    > qnorm(0.99, mean = 0, sd = 1)

Dando como resultado:

    [1] 2.326348

Se puede obtener más información sobre las distribuciones de probabilidad en el contexto de R con las funciones de ayuda: ?Normal, ?Binomial, ?TDist, ?Chi-squared, etcétera. Se sale de la página de ayuda con q.

Construcción de un Histograma

El archivo Histograma.Normal.r contiene los comandos necesarios para generar la siguiente figura.

En particular, el comando hist(x) genera el histograma de los datos $ n $ en el vector $ x $ .

   dev.new(width=6, height=5)
   par(cex=1.5)
   y <- rnorm(1500, 0, 10)
   hist(y, breaks=50, col="yellow", freq=FALSE)
   x <- seq(from=-40, to=40, by=1)
   w <- dnorm(x, mean=0, sd=10)
   lines(x, w, lwd=3, col="red")

Ley de los Grandes Números

En el histograma anterior, el tamaño de la muestra fue de 1500. Por la ley de los grandes números, cuando el tama ̃no de la muestra se incrementa, entonces el histograma tiende a coincidir con más exactitud con la función de densidad de la que se tomó la muestra.

Teorema del Límite Central

“As a rule of thumb, the sample size must be at least 30 for the central limit theorem to hold true”. Charles Wheela

“I know of scarcely anything so apt to impress the imagination as the wonderful form of cosmic order expressed by the law of frequency of error. The law would have been personified by the Greeks if they had known of it”. Francis Galton

Key Points

• First key point. Brief Answer to questions. (FIXME)

Distribuciones de Probabilidad que Dependenden de Uno o Más Parámetros

Overview

Teaching: 0 min
Exercises: 0 min

Questions

• Key question (FIXME)

Objectives

• First learning objective. (FIXME)

El científico realiza mediciones y las usa para encontrar fórmulas matemáticas que describan la naturaleza. Karl Pearson (1857-1936) tuvo la idea de que cuando se realiza un número grande de mediciones, lo que se obtiene es una distribución de valores. Esta distribución de valores se describe mediante una fórmula matemática que depende de uno o mas parametros.

La distribución Gama y sus parámetros

La función de densidad de una variable Gamma está dada por:

[\frac{1}{\Theta^{\alpha}\Gamma(\alpha)}x^{\alpha-1}e^{-x/\Theta}]

en donde $ \alpha $ es el parámetro de forma y $ \Theta $ el parámetro de escala. En la figura, $ \alpha $ toma los valores 1.5, 2, 3, 4.5, 6. El parámetro de escala $ \Theta $ se mantuvo constante e igual a 1. Ver archivo Gama.1.r.

Método del rechazo para generar muestras

El método del rechazo nos permite simular la realización de una variable aleatoria definida por una función de densidad arbitraria. Esta función puede depender de uno o más parámetros. Ver el archivo Método.del.Rechazo.r .

> ejectionK <- function(fx, a, b, K) {
> # simula una variable con funci ́on de densidad fx
> # supone que fx es 0 fuera de [1, b] y acatada por K
>   while (TRUE) {
>     x <- runif(1, a, b)
>     y <- runif(1, 0, K)
>     if (y < fx(x)) return(x)
>   }
> }

1. Genere X∼Uni[a,b] y Y∼Uni[0,K], variables independientes.
2. Si Y < fx(X) entonces se reporta X, en otro caso, ir al paso 1.

Método de Momentos

Los estimadores por el método de los momentos están dados por $ \widehat{\alpha}=\frac{m_{1}^{2}}{m_{2}-m_{1}^{2}} $ y $ \widehat{\Theta}=\frac{m_{2}-m_{1}^{2}}{m_{1}} $ en donde $ m_{k}= \frac{1}{n}\sum_{i=1}^{n}x_{i}^{k} $ .

En la figura se muestra el histograma de 10 000 estimaciones realizadas a partir de muestra de tamaño 50 cada una. La línea punteada en rojo muestra el valor real del parámetro.

Ver archivo Gama.Momentos.r .

Tamaño de la muestra:

> muestra <- 200

Número de muestras:

> tot <- 10000

Parámetro de forma:

> a <- 10

Parámetro de escala:

> t <- 0.5

Se inicializan las variables:

> t.hat <- c()
> a.hat <- c()

Y después:

> for(i in 1:tot) {
>   m <- rgamma(muestra, shape=a, scale=t)
>   m1 <- mean(m)
>   m2 <- sum(m^2)/muestra
>   t.hat[i] <- (m2-m1^2)/m1
>   a.hat[i] <- m1 ˆ 2/(m2-m1^2)
> }

Para ver la salida del código anterior:

> hist(a.hat)
> hist(t.hat)

La siguiente tabla contiene algunos estimadores según el método de los momentos. Note que $ m_{1}=\overline{X} $ .

La distribución logarítmica tiene distribución de masa y media dadas respectivamente por $ f(k)= \frac{-p^{k}}{klog(1-p)} $ y $ M_{1}(p)=\frac{-p}{(1-p)log(1-p)} $ .

Desafortunadamente no es posible resolver la ecuación $ M_{1}(p) = m_{1} $ para obtener $ \widehat{p} $ en función de $ m_{1} $ . Sin embargo, esta ecuación se puede resolver numéricamente. Para este fin, se define la siguiente función.

fun <- function(p) {
  a <- (1-p)∗log(1-p)
  return(-p/a-m1)
}

Con el siguiente comando se obtiene el valor numérico para $ \widehat{p} $ una vez que $ m_{2} $ es conocido.

uniroot(fun, c(0.01, 0.9))$root

Ver el archivo Distribución.Logarítmica.r. En este archivo también se programaron los comandos para realizar la simulación de la variable aleatoria.

> fd <- function(n=k, probabilidad=p) { # la distribución de masa
>  a <- -probabilidad^n
>  b <- n∗log(1-probabilidad)
>  return(a/b)
> }
> sim <- function(p) { # función que simula una variable logarítmica
>  rnd <- runif(1)
>  k <- 0
>  sum <- 0
>  while(rnd>sum) {
>     k <- k+1
>     sum <- sum+fd(k, p) 
>  }
> return(k)
> }

Se generan tot valores de la variable aleatoria:

tot <- 50000
u <- matrix(0.35, tot, 1)
z <- sapply(u, sim)

Se utilizan los valores generados para estimar el parámetro p:

print(uniroot(fun, c(0.01, 0.99), mean(z))$root)

Key Points

• First key point. Brief Answer to questions. (FIXME)

Métodos de Estimación de Parámetros, Puntual y por Intervalos

Overview

Teaching: 0 min
Exercises: 0 min

Questions

• Key question (FIXME)

Objectives

• First learning objective. (FIXME)

Máxima Verosimilitud

Ronald Fisher (1890-1962) se dio cuenta de que los métodos que Karl Pearson había estado usando para estimar los parámetros de una distribución, producían estadísticos que no necesariamente eran consistentes y que con frecuencia presentaban un sesgo. Para producir estadísticos consistentes y eficientes, Fisher propuso algo que ́él denominó «estimadores de máxima verosimilitud».

La función de verosimilitud de una muestra $ X_{1},…, X_{n} $ está dada por:

El estimador máximo verosímil $ \widehat{\Theta} $ del parámetro $ \Theta $ , es aquel valor para el cual:

[\mathscr{l}(\Theta)= log L(\Theta)]

es máximo. Si el tamaño $ n $ de la muestra es grande, entonces $ \widehat{\Theta} $ tiene una distribución aproximadamente normal con media $ \Theta $ y varianza:

Note que $ var(\widehat{\Theta_{n}})\to 0 $ cuando $ n\to \infty $ . Ver el archivo EMV.Exponencial.r.

> LL <- function(r) { #el log de la función de verosimilitud
>   R <- dexp(x, r)
>   -sum(log(R))
> }

> rt <- 2.5; N <- 1250; x <- rexp(N, rate=rt)
> optim(rt, LL, method = “Brent”, lower=0.001, upper=100)$par

> tot <- 2000
> a <- c()
> for(i in 1:tot) {
>     x <- rexp(N, rate=rt)
>     a[i] <- optim(rt, LL, method = "Brent", lower=0.001,
>                   upper=100)$par
>     if(i % %100==0) {print(i)}
> }

> hist(a)

Máxima Verosimilitud: Ajuste de una distribución

Con el comando fitdistr de la librería MASS es posible estimar los parámetros de una distribución mediante el método de máxima verosimilitud. He aquí algunos ejemplos.

> x <- rnorm(500, mean=0, sd=1)
> library(MASS)
> fitdistr(x, "normal")

> x <- rgamma(100, shape=10, scale=0.5)
> fitdistr(x, "gamma")

> x <- rgeom(100, p=0.1)
> fitdistr(x, "geometric")

Nota

En el último ejemplo, el comando fitdistr se ejecutó con el nombre “geometric” y no con el nombre “geom”. Consulte la siguiente página.

Extracción de información en objetos

Considere nuevamente la tarea del ajuste de una distribución

> x <- rnorm(500, mean=0, sd=1)
> library(MASS)
> fit <- fitdistr(x, "normal")

El comando summary(fit) despliega el contenido en el objeto fit.

Mediante fit$estimate, o fit[[1]], se obtienen las estimaciones obtenidas.

Mediante fit$estimate[[1]], o fit[[1]][[1]], se obtiene el valor estimado para el primero de los parámetros

Error Cuadrático Medio (ECM)

Un estimador $ \widehat{\Theta_{1}} $ es mejor que $ \widehat{\Theta_{2}} $ si $ ECM(\widehat{\Theta_{1}}) $ < $ ECM(\widehat{\Theta_{2}}) $ .

Si $ \widehat{\Theta} $ es insesgado, entonces $ ECM(\widehat{\Theta})= var(\widehat{\Theta}) $ .

Ejemplo. La población: $ X ∼ $ Uniforme $ (0,\Theta) $ . Los estimadores: $ \widehat{\Theta_{1}}= 2\overline{X} $ y $ \widehat{\Theta_{2}}= \frac{n +1}{n}max(X_{1},…,X_{n}) $ .

Las varianzas de los estimadores:

$ var(\widehat{\Theta_{2}})=\frac{\Theta^{2}}{n(n+2)} $ y $ var(\widehat{\Theta_{1}})=\frac{\Theta^{2}}{3n} $ .

Ver el archivo Contraejemplo.r.

El p-ésimo cuantil de una variable aleatoria $ X $ es aquel valor $ \phi_{p}=\phi_{p}(X) $ para el cual:

[p = P{ X \le \phi_{p}}]

Para calcular $ \phi_{0.75} $ de una Normal(0,1):

> qnorm(0.85, 0, 1)

[1] 1.036433

Intervalos de Confianza

Un intervalo $ (L,U) $ del $ (1-\alpha)\% $ de confianza para un parámetro $ \Theta $ es aquel para el cual:

[P{ L \le \Theta\le U}=1-\alpha]

Por ejemplo, el estimador de máxima verosimilitud para $ \Theta $ que se obtiene a partir de una muestra $ X_{1},…,X_{n} $ de la población $ X \sim $ Exponencial $ (\Theta) $ tiene una distribución apróximadamente normal

Por lo tanto, un intervalo del 0.95% de confianza está dado por:

[(\widehat{\Theta}+\frac{\widehat{\Theta}}{\sqrt{n}}\phi_{0.025},\widehat{\Theta}+\frac{\widehat{\Theta}}{\sqrt{n}}\phi_{0.975})]

con $ \phi_{p}=\phi_{p}(Z) $ .

En la figura se muestran 70 intervalos de confianza de nivel 95%.

Los intervalos fueron calculados usando la fórmula de la transparencia anterior. Algunos de los intervalos contienen al valor del parámetro verdadero (marcado en linea negra). Un porcentaje de aproximadamente 5 % de los intervalos, no contienen al parámetro verdadero. Ver archivo IntervaloDeConfianza.nb.

Key Points

• First key point. Brief Answer to questions. (FIXME)

El concepto de Prueba de Hipótesis Estadística y ejemplos de las pruebas más comunes.En particular las pruebas de bondad de ajuste.

Overview

Teaching: 0 min
Exercises: 0 min

Questions

• Key question (FIXME)

Objectives

• First learning objective. (FIXME)

Prueba de Hipótesis

Tres personajes son responsables de haber desarrollado el análogo estadístico al falsacionismo de Popper.

Una prueba de hipótesis es forma de verificar una afirmación sobre la distribución de una variable aleatoria $ X $ . Se trata de ver si la información acerca de $ X $ que está contenida en una muestra $ X_{1},…,X_{n} $ , tiende a confirmar la afirmación o más bien la contradice.

Consideremos un ejemplo. Sabemos que $ X \sim Normal(\mu,2) $ , en donde $ \mu\in{0,2 } $ . ¿Cómo podemos utilizar la información en una muestra $ X_{1},…,X_{n} $ , para decidir si la afirmación:

[H_{0}: \mu = 0]

es cierta o no?

Note que si $ H_{0} $ es cierta, entonces $ \overline{X} $ debería estar muy cerca de 0. ¿Qué pasa si observamos que $ \overline{X} = 2 $ ? Por otro lado, ¿es el valor $ \overline{X} = 0.75 $ evidencia suficiente como para afirmar que $ H_{0} $ es falsa?

Es natural rechazar la hipótesis $ H_{0} $ cuando $ \overline{X} $ sea grande, digamos, mayor que $ q $ . ¿Qué tan grande debe ser $ q $ para que la probabilidad

de equivocarnos cuando rechazamos $ H_{0} $ sea tolerable?

Cuando el tamaño de la muestra es igual a $ n $ , entonces $ \overline{X}\sim Normal(\mu,\frac{2}{\sqrt{n}}) $ . Poniendo $ \alpha = 0.05 $ se obtiene que

Por lo tanto $ q = \frac{2}{\sqrt{n}}\phi_{0.095}(Z)=\frac{3.28971}{\sqrt{n}}. $

El siguiente código implementa la prueba de hipótesis que se dedujo en la transparencia anterior

> mu <- 0  #el valor de la media de la población
> n <- 5  #el tamaño de la muestra
> q <- 3.289707/sqrt(n) #ver ejercio anterior
> tot <- 10  #número de veces que se aplica la prueba
> sum <- 0 
> for(i in 1:tot) {
>   x <- rnorm(n, mean=mu, sd=2)
>   x.bar <- mean(x)
>   if(x.bar < q) {sum <- sum+1}
> }
> print(sum/tot) #porcentaje de error tipo I

Córrase este código con distintos valores para la variable tot. Ver el archivo Prueba.de.Hipótesis.r .

Prueba de Hipótesis: Tipos de Error

Error tipo I: rechazar $ H_{0} $ cuando $ H_{0} $ correcta.

Error tipo II: no rechazar $ H_{0} $ cuandob$ H_{0} $ incorrecta.

El falsacionismo de Karl Popper

El falsacionismo, o principio de falsabilidad, es una corriente epistemológica fundada por Karl Popper (1902-1994). Para Popper, contrastar una teoría significa intentar refutarla mediante un contraejemplo. Si no es posible refutarla, dicha teoría queda corroborada, pudiendo ser aceptada provisionalmente, pero no verificada; es decir, ninguna teoría es absolutamente verdadera, sino a lo sumo se mantiene como no refutada. El falsacionismo es uno de los pilares del método científico.

Prueba de hipótesis para la diferencia de medias.

Suponga que tenemos dos poblaciones:

[X \sim Normal(\mu_{x},\sigma_{x}^{2}) \; y \; Y \sim Normal(\mu_{y},\sigma_{y}^{2})]

Estamos interesados en probar la hipótesis:

[H_{0}: \mu_{x}-\mu_{y}=\Delta_{0} \; conta \; H_{1}: \mu_{x}-\mu_{y}\neq \Delta_{0}]

Sean $ X_{1},…,X_{n_{x}} $ y $ Y_{1},…,Y_{n_{y}} $ muestras de las poblaciones $ X $ y $ Y $ . Es estadístico de prueba es:

$ T_{0}=\frac{\overline{X}-\overline{Y}-\Delta_{0}}{\sqrt[s_{p}]{\frac{1}{n_{x}}+\frac{1}{n_{y}}}} $ en donde $ s_{p}^{2}= \frac{(n_{x}-1) + (n_{y}-1)s_{y}^{2}}{n_{x}-n_{y}-2} $ y $ s_{x}^{2}= \frac{1}{n_{x}-1}\sum_{i=1}^{n_{x}}(X_{i}-\overline{X})^{2} $ el estimador de la varianza $ \sigma_{x}^{2} $ .

Es claro que valores de $ \abs{T_{0}} $ pequeños son congruentes con la hipótesis $ H_{0} $ . ¿Qué tan grande debe ser $ \abs{T_{0}} $ para poder rechazar $ H_{0} $ ? Puesto que la distribución del estadístico de prueba es conocida (distribución T de Student) entonces es posible calcular el valor $ q $ tal que $ H_{0} $ se rechaza siempre que $ \abs{T_{0}} \gt q $ con una probabilidad de error tipo I prefijada.

El comando:

> t.test(x, y, mu=0.0)

Realiza la prueba de la hipótesis $ H_{0} $ en donde $ mu = \Delta_{0} $ . Como resultado de ejecutar este comando se obtendrá el así llamado p-valor.

William Sealy Gosset

William Sealy Gosset (1876-1937) fue un estadístico, mejor conocido por su sobrenombre literario Student. Estudió química y matemática en el New College de Oxford. Tras graduarse en 1899, se incorporó a las destilerías Guinness en Dublín. Guinness era un negocio agroquímico y ahí Gosset pudo aplicar sus conocimientos estadísticos tanto a la destilería como a la granja para seleccionar las mejores variedades de cebada.

Gosset introdujo la distribución T de Student para realizar la prueba de diferencia de medias.

[T(k)=\frac{Z}{\frac{Ji^{2}(k)}{k}}]

El p-valor de una prueba de hipótesis

El p-valor es la probabilidad de que se observe un valor del estadístico de prueba como el que se ha observado al realizar el experimento, bajo el supuesto de que $ H_{0} $ es cierta.

Por lo tanto, p-valores pequeños deben interpretarse como evidencia en contra de $ H_{0} $ .

Como ayuda para asimilar el concepto de p-valor, se puede correr el código en el archivo Prueba.T.de.Student.r con distintos para las medias de $ X $ y $ Y $ .

Tamaño de la muestra de $ X $ :

> nx <- 100

Tamaño de la muestra de $ Y $ :

> ny <- 50

Desviación estándar común de las dos poblaciones:

> desviacion <- 1
> x <- rnorm(nx, 0.0, sigma)
> y <- rnorm(ny, 0.5, sigma)

Prueba de la hipótesis $ H_{0}: \mu_{x}= \mu_{y} $ :

> t.test(x, y, mu=0)

Aplicación al consumo de energía eléctrica

En la figura se muestran dos histogramas sobrepuestos. El histograma en color rojo representa el consumo personal diario de energía eléctrica durante los meses de verano, y en verde, en consumo en los meses de invierno.

Se aplicó la prueba T de Student para contrastar la hipótesis nula según cual no existe diferencia en los consumos, y se observó un p-valor igual a 0.02701, con lo cual podemos rechazar esta hipótesis. Ver el archivo Consumo.Electricidad.r.

Prueba para la igualdad de dos varianzas

Sean $ \sigma_{1}^{2} $ y $ \sigma_{2}^{2} $ las varianzas de dos poblaciones normales e independientes. Para probar la hipótesis:

[H_{0}:\sigma_{1}^{2}=\sigma_{2}^{2}]

se usa el comando var.test cuyos dos primeros argumentos son dos vectores numéricos que contienen los datos de cada muestra.

> x <- rnorm(50, mean=0, sd=1)
> y <- rnorm(50, mean=0, sd=2)
> var.test(x, y)

F test to compare two variances data: x and y

$ F = 0.28962 $ , $ num df = 49 $, $ denom df = 49 $ , $ p-value = 2.808e-05 $

alternative hypothesis: true ratio of variances is not equal to 1

Nota

Ver la página siguiente para más detalles de este comando.

var.test: F Test to Compare Two Variances.

Prueba de significancia para la correlación

El comando cor.test(x, y) nos permite probar la hipótesis nula de que los vectores $ x, y $ no están correlacionados.

Definimos la relaci ́on entre $ x, y $ :

> recta <- function(x) 1+5∗x 
> set.seed(415)
> x <- c(1:10)
> y <- recta(x)+rnorm(10, 0, 5)
> cor.test(x, y, method="pearson") #otro método: spearman

Pearson’s product-moment correlation data: x and y

t = 6.0949, df = 8, p-value = 0.0002911

alternative hypothesis: true correlation is not equal to 0

95 percent confidence interval: 0.6469479 0.9780966

sample estimates: cor 0.907086

La prueba de bondad de ajuste de Kolmogorov-Smirnov

Dada una muestra $ X_{1},…,X_{n} $ , tomada de una población $ X $ , nos interesa probar la siguiente hipótesis:

[H_{0}: X\;tiene\;una\;distribución\;F(x)]

El comando ks.test nos permite realizar la prueba de bondad de ajuste de Kolmogorov y Smirnov. He aquí algunos ejemplos.

> y <- rnorm(20, mean=0, sd=5)
> ks.test(y, "pnorm", mean=0, sd=3)

> y <- rchisq(50, df=5)
> ks.test(y, "pchisq", df=2)

Si el p-valor obtenido es pequeño, entonces no es creíble que la muestra haya sido generada por la distribución $ F(x) $ .

La transformación Box-Cox

Cuando un conjunto de datos no se distribuye normalmente, entonces es posible transformarlos de manera que los datos transformados sí se distribuyan normalmente. En el lado derecho tenemos una muestra de una población gamma. En el lado derecho tenemos la misma muestra despu ́es de que se aplicó la transformación Box-Cox. Ver el archivo Box.Cox.r.

Aquí se puede usar la prueba de Kolmogorov-Smirnov para evaluar el ajuste de los datos transformados a la distribución normal.

La prueba Ji cuadrada de bondad de ajuste

Suponga que una variable $ X $ tiene una función de distribución de masa de probabilidad (dmp) dada por la siguiente tabla.

Dada una muestra de la cual sospechamos que fue tomada de la distribución dmp, con el comando:

> chisq.test(muestra, p=dmp)

Se obtiene el p valor para la prueba de la hipótesis nula $ H_{0} $ que afirma que la muestra fue tomada de la población $ X $ .

Ver el archivo Ji.Cuadrada.r.

La prueba $ X^{2} $ de bondad de ajuste

En el siglo XIX, se aplicaron m ́etodos estad ́ısticos a datos biológicos. Los investigadores supon ́ıan que los datos seguían una distribución normal. Pero, en 1900, Karl Pearson criticó este supuesto y observó que los histogramas obtenidos presentaban una asimetría (skewness).

En una serie de artículos entre 1893 y 1916, Pearson introdujo una familia de distribuciones, de la cual, la distribución normal era un caso particular. Pearson introdujo la prueba Ji cuadrada para mostrar que las muestras obtenidas no no podían provenir de una distribución normal.

Key Points

• First key point. Brief Answer to questions. (FIXME)

Regresión Lineal

Overview

Teaching: 0 min
Exercises: 0 min

Questions

• Key question (FIXME)

Objectives

• First learning objective. (FIXME)

Regresión Lineal Simple

Considere una variable $ Y $ de la forma

[Y = \beta_{0} + \beta_{1} x + \epsilon]

en donde $ \beta_{0} $ y $ \beta_{1} $ son dos parámetros del modelo, $ x $ es una variable no aleatoria y

[\epsilon \sim Normal(0, \sigma^{2})]

Aquí, $ \sigma^{2} $ es el tercer parámetro del modelo. Suponga que tenemos n pares de datos

[(x_{1}, Y_{1}), (x_{2}, Y_{2}), \cdots, (x_{n}, Y_{n})]

A partir de esta muestra, se obtienen estimaciones para los tres parámetros del modelo $ \hat{\beta}_{0} $, $ \hat{\beta}_{1} $ y $ \hat{\sigma}^{2} $

Ejemplo de aplicación del modelo

Propiedades de los estimadores del modelo

Notación

$ e_{i} = y_{i} - \hat{y}_{i} $ es el i-ésimo residual, siendo $ \hat{y}_{i} = \hat{\beta}_{0} + \hat{\beta}_{1} x_{i} $

$ SCE = \sum\limits_{i=1}^{n} (y_{i} - \hat{y}_{i})^{2} $ Es la suma de los cuadrados del error

1.     $ \hat{\beta}_{0} = \bar{y} - \hat{\beta}_{1} \bar{x} $        y        $\hat{\beta}_{1} = \frac{S_{xy}}{S_{xx}} $        en donde        $ S_{xy} = \sum\limits_{i=1}^{n} (x_{i} - \bar{x}) (y_{i} - \bar{y}) $
   
   
2.     $ E(\hat{\beta}_{0}) = \beta_{0} $        y        $ E(\hat{\beta}_{1}) = \beta_{1} $
   
   
3.     $ \textrm{var}(\hat{\beta}_{1}) = \frac{\sigma^{2}}{S_{xx}} $        $ \textrm{var}(\hat{\beta}_{0}) = (\frac{1}{n} + \frac{\bar{x}^{2}}{S_{xx}}) \sigma^{2} $
   
   
4.     $ \hat{\beta}_{1} $,     $ \bar{Y} $     y     $ SCE $     son independientes
   
   
5.     $ \sum\limits_{i = 1}^{n} (\frac{e_{i}}{\sigma})^{2} \sim Ji^{2}(n-2) $        y        $ \hat{\sigma}^{2} = \frac{SCE}{n-2} =: $ CME es insesgado

Las distribuciones de la teoría estadística

\(Z = \frac{1}{\sqrt{\textrm{var}(x)}} \sum\limits_{i = 1}^{n} (X_{i} - \mu)\)     (Teorema del límite central)

\(Ji^{2}(k) = \sum\limits_{i = 1}^{k} Z_{i}^{2}\)     (Suma de cuadrados normales estándar)

\(T(k) = \frac{Z}{\sqrt{Ji^{2} (k) / k}}\)     (Estandarización usando $\hat{\sigma}$ en lugar de $\sigma$)

\(F(k, r) = \frac{Ji^{2}(k)/k}{Ji^{2}(r)/r}\)     (Cociente de Ji cuadradas)

Francis Galton

Francis Galton (1822-1911) fue un polímata, antropólogo, geógrafo, explorador, inventor, meteorólogo, estadístico, psicólogo y eugenista británico. Creó el concepto estadístico de correlación y regresión hacia la media, altamente promovido. El fue el primero en aplicar métodos estadísticos para el estudio de las diferencias humanas y la herencia de la inteligencia, introdujo el uso de cuestionarios y encuestas para recoger datos sobre las comunidades humanas.

Regresión hacia la media

En estadística, la regresión hacia la media es el fenómeno en el que si una variable es extrema en su primera medición, tenderá a estar más cerca de la media en su segunda medición y, paradójicamente, si es extrema en su segunda medición, tenderá a haber estado más cerca de la media en su primera.

Regresión lineal simple

Diagrama de dispersión

recta <- function(x) 1 + 5*x  #Defina una función lineal "recta(x)"
set.seed(415)
x <- c(1:10)
y <- recta(x) + rnorm(10, 0, 5)
dev.new(width=6, height=5)
par(cex=1.5)
plot(x, y)  #Para obtener el diagrama de dispresión

Recta estimada

cor(x, y) # calcula la correlación entre x, y

[1] 0.907086

A <- lm(y ∼ x) # asigna a A el modelo estimado y=b+m∗x
print(A) # despliega los coeficientes calculados para el modelo

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
      5.192        4.666 

abline(A) # añade la recta de regresión al diagrama de dispersión

plot(A, which=1, add.smooth=FALSE) # grafica los residuales

Estadísticos del modelo

summary(A) # para obtener el resumen del modelo

Call:
lm(formula = y ~ x)

Residuals:
   Min     1Q Median     3Q    Max 
-7.623 -3.364 -1.708  3.693 14.355 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.1919     4.7498   1.093 0.306184    
x             4.6657     0.7655   6.095 0.000291 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.953 on 8 degrees of freedom
Multiple R-squared:  0.8228,	Adjusted R-squared:  0.8007 
F-statistic: 37.15 on 1 and 8 DF,  p-value: 0.0002911

Coeficiente de determinación $ R^{2} $

El coeficiente de determinación de la regresión está dado por:

[R^{2} = \frac{\textrm{suma de cuadrados de la regresión}}{\textrm{suma de cuadrados del error}} = \frac{S_{xy}^{2}}{S_{xx} S_{yy}}]

El estadístico $R^{2}$ nos dice cuál es el porcentaje de variación de los datos que se explica por el modelo. $R^{2}$ no debe usarse para evaluar el ajuste del modelo a los datos

Gráfica de los residuales

set.seed(405)
x <- seq(from=0.1, to=10, by=0.1)
Y <- recta(x)+rnorm(length(x), 0, 5)
A <- lm(Y ∼ x)
modelo <- function(x) -0.3041+5.2817∗x
Y1 <- modelo(x)
qqnorm(Y-Y1, pch="o")

qqline(Y-Y1)

Estructura de datos

R usa cuatro tipos básicos de estructuras de datos:

1. Un vector puede tener entradas numéricas, de texto o de valores lógicos.
2. Un arreglo (array) es un vector con especificaciones de dimensión. El arreglo más común es la matriz.
3. Un factor es un objeto usado para definir variables categóricas.
4. Un marco de datos (data frame) es una lista de vectores de la misma dimensión. Cada vector en un marco de datos corresponde a una variable de un experimento

vec2 <- c(1, 2, 3.4, 5.6, 7)
vec3 <- c("tipo A", "tipo B", "tipo C", "tipo D", "tipo E")
vec4 <- c(TRUE, TRUE, FALSE, TRUE, FALSE)
f1 <- factor(vec3) # usa el vector vec3 para hacer el factor f1
datos <- data.frame(vec2, f1, vec4) 
print(datos)

  vec2     f1  vec4
1  1.0 tipo A  TRUE
2  2.0 tipo B  TRUE
3  3.4 tipo C FALSE
4  5.6 tipo D  TRUE
5  7.0 tipo E FALSE

Con el comando fix(datos) se pueden hacer modificaciones a este marco de datos

Regresión Lineal Multiple

set.seed(777)
f <- function(x, y) 3∗x-5∗y+1
x <- runif(20, 0, 10)
y <- runif(20, 0, 10)
resp <- f(x, y) + rnorm(20, 0, 3)
datos <- data.frame(x, y, resp)
plot(datos, cex=1.5, cex.axis=2) # genera la gr´afica: >lm(datos$resp∼datos$x+datos$y)

lm(datos$resp∼datos$x+datos$y)

Call:
lm(formula = datos$resp ~ datos$x + datos$y)

Coefficients:
(Intercept)      datos$x      datos$y  
     -0.658        3.256       -4.958 

Ver el archivo Regresión.Múltiple.r

Varios comandos para la regresión lineal

Suponga que se tienen datos de la siguiente forma

[y_{i} = \beta_{0} + \beta_{1} u_{i} + \beta_{2} v_{i} + \beta_{3} w_{i} + \epsilon_{i}]

Con el comando >fit <- lm(y ∼ u + v + w) se obtiene el ajuste del modelo. La tabla siguiente contiene comandos para obtener distintas estadísticos de la regresión

Ver el archivo Regresión.Comandos.Varios.r

Con el comando >lm(y ∼ x + 0) se ajusta el modelo sin término constante $ y_{i} = \beta x_{i} + \epsilon_{i} $

Con el comando >lm(y ∼ u∗v) se ajusta el modelo

[y_{i} = \beta_{0} + \beta_{1} u_{i} + \beta_{2} v_{i} + \beta_{3} u_{i} v_{i} + \epsilon_{i}]

Ver el archivo Regresión.Términos.Cruzados.r

Con los dos comandos que siguen, es posible seleccionar el mejor conjunto de variables explicativas.

full.model <- lm(y ∼ u + v + w + z)
reduced.model <- step(full.model, direction="backward")

Ver el archivo Regresión.Modelo.Reducido.r

Análisis de regresión

Variables indicadoras

Considere el modelo $ y = \beta_{0} + \beta_{1} u+ \beta_{2} v + \epsilon $, en donde la variable $ v $ solamente puede tomar los valores 0 o 1. En la gráfica los valores de $ y $ que se obtienen cuando $ v = 0 $, aparecen en negro y los mismos valores cuando $ v = 1 $, aparecen en rojo. Es interesante la prueba de hipótesis según la cual dos interceptan (cuando $ v = 0 $ y $ v = 1 $) son iguales

Ver el archivo Variables.Indicadores.r

Con el modelo de la lámina anterior, $ Y = \beta_{0} + \beta_{1} u + \beta_{2} v + \epsilon $ consideramos ahora el ajuste del modelo ampliado

[Y = \beta_{0} + \beta_{1} u + \beta_{2} v + \beta_{12} u v + \epsilon]

Cuando $ v = 0 $, se obtiene que

[E(Y) = \beta_{0} + \beta_{1} u]

Cuando $ v = 1 $, se obtiene que

[E(Y) = \beta_{0} + \beta_{1} u + \beta_{2} + \beta_{12} u = (\beta_{0} + \beta_{2}) + (\beta_{1} + \beta_{12}) u]

Son interesante las hipótesis nulas $ H_{0} : \beta_{2} = 0 $ y $ H_{0} : \beta_{12} = 0 $

Regresión polinomial

Cuando queremos ajustar a un conjunto de datos, un modelo de la forma

[y = \beta_{0} + \beta_{1} u + \beta_{2} u^{2} + \cdots + \beta_{p} u^{p} + \epsilon]

se usa el comando >lm(y ∼ poly(u, p, raw=TRUE)) , en donde $ p $ es el grado del polinomio que queresmo ajustar.

Ver el archivo Regresión.Polinomial.r

Key Points

• First key point. Brief Answer to questions. (FIXME)

Análisis de Varianza

Overview

Teaching: 0 min
Exercises: 0 min

Questions

• Key question (FIXME)

Objectives

• First learning objective. (FIXME)

Análisis de Varianza

En 1919, la Estación Rothamsted contrató al joven estadístico Ronals Aylmer Fisher para que aprovecara los datos ahí acumulados. El análisis de Fisher sugería que la relación entre la lluvia y el crecimiento de las plantas era más significativa que la relación entre el fertilizante y el crecimiento de las plantas. A los científicos de la estación de Rothamsted, no les interesaba la lluvia como factor determinante de la cosecha, sino el fertilizante. No se sabía separar los efectos d ela lluvia d elso efectos del fertilizante. Fisher comprendió que los efectos se podían seprar si los experimentos se diseñaban de manera apropiada.

El análisis de varianza (ANOVA) de un sólo factor se utiliza para comparar las medias I poblaciones. Es de interés la hipótesis

[H_{0} : \mu_{1} = \mu_{2} = \cdots = \mu_{I}]

En el ANOVA se tienen datos de la forma

[X_{ij} = \mu_{i} + \epsilon_{ij}]

en donde $\epsilon_{ij}$ son las desviaciones (o errores) aleatorios que la j-ésima observación respecto de la i-ésima media poblacional $\mu_{i}$. Se supone que los errores $\epsilon_{ij}$ son independientes con media cero y desviación estándar constante.

Si $ H_{0} $ se cumple, entonces el estadístico de prueba $ F_{cal} = CMTr/CME $ tiene un valor de 1. Valores grandes de $ F_{cal} $ sugieren que $ H_{0} $ es falsa. ¿Qué tan grande debe ser $ F_{cal} $ para poder rechazar a $ H_{0} $?

set.seed(777)
a <- 1.0 + rnorm(10, 0, 0.7)
b <- 1.5 + rnorm(10, 0, 0.7)
c <- 2.0 + rnorm(10, 0, 0.7)
obs <- c(a, b, c)
A <- rep("A", 10)
B <- rep("B", 10)
C <- rep("C", 10)
fac <- factor(c(A, B, C))
datos <- data.frame(obs, fac)
mod <- oneway.test(obs ∼ fac, data=datos, var.equal=TRUE)
print(mod) # se obtiene como resultado lo siguiente:

	One-way analysis of means

data:  obs and fac
F = 1.6901, num df = 2, denom df =
27, p-value = 0.2034

En la figura de la izquierda se reportan los diagramas de caja de las tres poblacionse generadas por el código anterior. Dentro de cada caja yacen las observaciones entre loso cuartiles $ Q_{1} $ y $ Q_{3} $ de cada población. Los “bigotes” se extienden hasta los valores máximo de la serie o hasta $ 1.5 x (Q_{3} - Q{1}) $.

En la figura de la derecha se reportan los intervalos “estudentizados” mediante el procedimiento de Tukey.

En estos diagramas de caja se muestra el consumo de personal de energía eléctrina por bimestre.

tabla

Para la hipótesis nula según la cual no hay diferencias en el consumo, se tiene un p-valor de 0.0126. Por lo tanto se rechaza esta hipótesis. En la figura, los intervalos estudentizados de Tukey.

Se observa la mayor diferencia entre los bimestres junio - julio contra febrero - marzo.

nota

El ANOVA de un factor como un caso de la regresión

En el archivo Solo.Variables.Indicadoras se considera un modelo de regresión de la forma

[Y = \beta_{0} + \beta_{1} v + \beta_{2} w + \epsilon]

en donde $ v $ y $ w $ son variables indicadoras que solamente toman los valores 0 a 1, pero no se da el caso de que $ v = w = 1 $. Note que

[E(Y) = \begin{cases} \beta_{0} &\quad\text{si } v = 0 \text{ y } w = 0
\beta_{0} + \beta_{1} &\quad\text{si } v = 1 \text{ y } w = 0
\beta_{0} + \beta_{2} &\quad\text{si } v = 0 \text{ y } w = 1
\end{cases}]

Cuando se varían los valores $ \beta_{0} $, $ \beta_{1} $, $ \beta_{2} $ se cambian los niveles del factor.

ANOVA con dos factores

Un acumulador de energía eléctrica se puede construir de tres tipo distintos de materiales (factor 1 con tres niveles) y trabajará a distintas temperaturas (factor 2). Es de interés el número de horas de servicio del acumulador. El acumulador será probado a tres temperaturas distintas (los tres niveles del factor 2). En la tabla se reportan las horas de servicio obtenidas al aplicar cada tratamiento (combinación de niveles de los factores) a 4 acumuladores elegidos de forma aleatoria

Con los siguientes comandos se realiza el análisis de varianza con dos factores para el problema planteado anteriormente.

vida <- c(130, 74, 150, 159, 138, 168, 155, 180, 188, 126, 110, 160, 34, 80, 135, 106, 174, 150, 40, 75, 122, 115, 120, 139, 20, 82, 25, 58, 96, 82, 70, 58, 70, 45, 104, 60)
material <- rep(c(1,1,2,2,3,3), 6)
temp <- rep(c(15, 70, 125), each=12)
dat <- data.frame(resp=vida, tipo=factor(material), temperatura=factor(temp))
fit <- aov(resp ∼ tipo∗temperatura, data=dat)
summary(fit)

            Df Sum Sq Mean Sq F value   Pr(>F)    
tipo         2  10678    5339   6.269  0.00531 ** 
temp         1  39043   39043  45.841 1.66e-07 ***
tipo:temp    2   2315    1158   1.359  0.27226    
Residuals   30  25551     852                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ver el archivo ANOVA.Acumulador.r

Regresión logística

Con los siguientes comandos se simula una muestra en donde la variable repuesta solamente toma los valores 0 o 1.

tt <- 50 # tama˜no de la muestra
x <- 20+c(1:tt)∗50/tt
b0 <- -10
b1 <- 0.2
y <- c()
for (i in 1:tt) {
  EXP <- exp(b0+b1∗x[i])
  p <- EXP/(1 + EXP)
  y[i] <- rbinom(1, size=1, prob=p)
}
plot(x, y)

Ver el archivo Regresión.Logística.r

Una vez dada la muestra que se generó en al página anterior, los coeficientes del modelo se estiman con el siguiente comando

fit <- glm(y ∼ x, family=“binomial”)
summary(fit)

En la figura se reporta la muestra generada junto con el modelo ajustado (en linea roja punteada). El modelo provee la probabilidad de que un punto de la muestra haya tomado el valor 1. Cuando esta probabilidad es muy pqeueña, entonces la variable respuesta tenderá a asumir el valor 0.

Key Points

• First key point. Brief Answer to questions. (FIXME)