Modelos Estadísticos 1

Obteniendo resultados estadísticos con R, parte IV

Posted in Tutoriales, Uncategorized by hector on mayo 19, 2009

Esta es la cuarta parte de una breve introducción a los comandos del sistema estadístico R. En esta parte trabajaremos con algunas ideas fundamentales del análisis estadístico y de prueba de hipótesis.


Regresión y Correlación

En muchas ocasiones tendremos variables (numéricas) que están relacionadas de manera lineal. Además, es importante notar que si las variables no tienen una relación lineal, dicha relación puede asumirse en una región de interés (i.e. se toma una aproximación lineal). El modelo de regresión lineal está dado por:
\displaystyle y_i = \alpha + \beta x_i + \varepsilon_i
donde 1 \leq i \leq N. Noten que x_i, y_i representan los distintos valores de las variables x,y; \varepsilon_i representa los errores asociados con dichas mediciones. Estos errores se asumen independientes y obtenidos de una distribución normal N(0,\sigma^2). Asociando esta expresión con la ecuación punto-pendiente de una recta y = m x + b, tenemos que \beta es la pendiente y se denota coeficiente de regresión; además \alpha es el intersecto. Nota: los parámetros \alpha, \beta, \sigma pueden estimarse empleando el Método de Mínimos Cuadrados en donde el objetivo es minimizar el cuadrado de los residuos:
\displaystyle E = \sum_{i=1}^{N} \left(y_i - (\alpha + \beta x_i)\right)^2
después de un poco de Cálculo llegamos a las expresiones:
\hat \beta = \frac{\sum (x_i - \bar x)(y_i - \bar y)}{\sum (x_i - \bar x)^2}
\hat \alpha = \bar y - \hat \beta \bar x

Un concepto bastante relacionado a la idea de regresión lineal es el de correlación lineal. Un coeficiente de correlación es una medida de la asociación entre dos variables; puede tomar valores entre -1 y +1 (donde 0 significa que no hay correlación entre las variables). Nota: “Correlación no implica causalidad”. La expresión para el coeficiente de correlación de Pearson (o r de Pearson) es:
r = \frac{\sum (x_i - \bar x)(y_i - \bar y)}{\sqrt{\sum(x_i - \bar x)^2 \sum(y_i - \bar y)^2}}

Veamos ahora algunos ejemplos en R. El comando para obtener la r de Pearson (coeficiente de correlación) es cor(x,y):

> x = seq(0, 10, by=0.1)
> y = 5*x + 7
> plot(x,y)
> cor(x,y)
[1] 1
> y = -5*x + 7
> plot(x,y)
> cor(x,y)
[1] -1

Agreguemos un poco de ruido (incertidumbre) a este ejemplo:

> x = seq(0, 10, by=0.1)
> y = 5*x + 7
> ruido = rnorm(length(x), mean=0, sd=3)
> plot(x, y + ruido)
> cor(x, y + ruido)
[1] 0.9818713

Vemos que el coeficiente de correlación (r de Pearson, en este caso) ya no es exactamente 1 debido al ruido introducido; sin embargo, aún existe una alta asociación (0.9818) entre las variables.

Los comandos importantes para regresión lineal son lm (linear model), abline. Se ejemplifican a continuación:

> x = seq(0, 10, by=0.1)
> y = 5*x + 7
> ruido = rnorm(length(x), mean=0, sd=3)
> plot(x, y + ruido)
> lm(y~x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
          7            5  

> lm(y + ruido ~ x)

Call:
lm(formula = y + ruido ~ x)

Coefficients:
(Intercept)            x  
      7.498        4.913  

> yr = y + ruido
> summary(lm(yr ~ x))

Call:
lm(formula = yr ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.0060 -1.7720 -0.4454  1.8719  6.3629 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.4980     0.5383   13.93   <2e-16 ***
x             4.9134     0.0930   52.83   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 2.725 on 99 degrees of freedom
Multiple R-squared: 0.9657,     Adjusted R-squared: 0.9654 
F-statistic:  2791 on 1 and 99 DF,  p-value: < 2.2e-16 

> plot(x,yr)
> abline(lm(yr ~ x))
> plot(x,yr)
> abline(7.4980, 4.9134)

regr-lin


Prueba t de una muestra

Empecemos dejando claro que la prueba t está diseñada en el supuesto que los datos vienen de una distribución normal N(\mu, \sigma^2). Podemos entonces estimar los parámetros \mu, \sigma de la distribución muestral a través de la media \overline{x} y desviación estándar s. La prueba t se emplea para probar hipótesis acerca de las medias de dos muestras, e.g. para probar si existe alguna diferencia entre el grupo control y el grupo experimental. Veamos algunos ejemplos.

Dieta recomendada.

Unos investigadores le han pedido a 11 estudiantes de secundaria que registren sus hábitos alimenticios. En base a las observaciones de los estudiantes se tienen los siguientes datos de consumo energético promedio en kilo-joules (1 kJ = 239 Cal). La dieta recomendada tiene un valor energético de 7725 kJ. ¿Qué tanto se desvían las dietas de las estudiantes de la dieta recomendada?

> consumo = c(5260,5470,5640,6180,6390,6515,6805,7515,7515,8230,8770)
> summary(consumo)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   5260    5910    6515    6754    7515    8770 
> t.test(consumo, mu=7725)

        One Sample t-test

data:  consumo 
t = -2.8208, df = 10, p-value = 0.01814
alternative hypothesis: true mean is not equal to 7725 
95 percent confidence interval:
 5986.348 7520.925 
sample estimates:
mean of x 
 6753.636 

La primera línea que contiene información que necesitamos interpretar es:

t = -2.8208, df = 10, p-value = 0.01814

En alguna época hubiéramos tenido que calcular t y luego buscar en una tabla de una distribución t el valor de p. Acá podemos leerlo inmediatamente. Dado que p < 0.05, tenemos entonces que rechazar la hipótesis nula: la dieta promedio de los estudiantes es igual a la dieta recomendada (ojo a la línea que despliega la hipótesis alternativa) Nota: un intervalo de confianza siempre se calcula seleccionando primero un nivel de confianza (parámetro conf.level por defecto es 0.95).


Consumo de gasolina.

La compañía de carros X ha anunciado que la nueva SUV tiene un consumo de gasolina de 17 millas por galón. Los editores de una revista automovilística consideran que la compañía X ha mentido y que dicha SUV hace menos millas por galón. Toman las siguientes mediciones y corren una prueba t:

> mpg = c(11.4,13.1,14.7,14.7,15.0,15.5,15.6,15.9,16.0,16.8)
> t.test(mpg, mu=17, alt='less')

        One Sample t-test

data:  mpg 
t = -4.2847, df = 9, p-value = 0.001018
alternative hypothesis: true mean is less than 17 
95 percent confidence interval:
     -Inf 15.78127 
sample estimates:
mean of x 
    14.87 

Dado que p-value < 0.05, nuevamente rechazamos la hipótesis nula ("la nueva SUV hace 17 millas por galón") en favor de la hipótesis alternativa ("la nueva SUV hace menos de 17 millas por galón"). Noten el parámetro alt='less' que modifica el análisis y la hipótesis alternativa.


Costos de cafetería.

La cafetería de la UVG dice que, en promedio, los estudiantes gastan Q.101.75 semanales en refacciones, almuerzos, etc. Un grupo de estudiantes de estadística realiza un muestreo aleatorio y encuentran los siguientes consumos semanales: 140, 125, 150, 124, 143, 170, 125, 94, 127, 53. La hipótesis nula en este caso es que el consumo semanal por estudiante es de Q101.75. La hipótesis alternativa es que el consumo semanal es mayor a lo que dice la cafetería.

Nota: dado que la prueba t asume que los datos vienen de una distribución normal, lo primero que deberíamos hacer es visualizar qué tan normales son los datos. Pueden emplear el comando qqnorm (vean su documentación). Un ejemplo:

> qqnorm(rnorm(1000))

Ahora la corrida de la prueba t:

> gasto = c(140, 125, 150, 124, 143, 170, 125, 94, 127, 53)
> qqnorm(gasto)
> t.test(gasto, mu=101.75, alt='greater')

        One Sample t-test

data:  gasto 
t = 2.291, df = 9, p-value = 0.02385
alternative hypothesis: true mean is greater than 101.75 
95 percent confidence interval:
 106.4169      Inf 
sample estimates:
mean of x 
    125.1 

Dado el valor p obtenido (i.e., p-value < 0.05), nuevamente rechazamos la hipótesis nula; en otras palabras: en base a los datos recopilados, el consumo semanal (por estudiante) en la cafetería es mayor a Q101.75. Nuevamente noten el parámetro alt='greater' y cómo modifica la hipótesis alternativa.


Análisis de Varianza

Supongamos que tenemos varias mediciones separados en distintos grupos. Empecemos con un poco de notación: x_{i,j} representa entonces la observación (o medición) j en el grupo i. Ahora un poco de manipulación algebraica:
\displaystyle x_{i,j} = x_{i,j}
\displaystyle = \bar x_i + x_{i,j} - \bar x_i
\displaystyle = \bar x + \bar x_i - \bar x + x_{i,j} - \bar x_i
\displaystyle = \bar x + (\bar x_i - \bar x) + (x_{i,j} - \bar x_i)

donde (\bar x_i - \bar x) es la desviación de la media del grupo (\bar x_i) con respecto a la media global (\bar x); (x_{i,j} - \bar x_i) es la desviación de la observación con respecto a la media del grupo. Esto corresponde entonces al modelo estadístico:
\displaystyle X_{i,j} = \mu + \alpha_i + \varepsilon_{i,j}
donde \varepsilon_{i,j} \sim N(0, \sigma^2). Tenemos entonces las siguientes definiciones:

  • Variación dentro de grupos: \displaystyle SS_W = \sum_i \sum_j (x_{i,j} - \bar x_i)^2
  • Variación entre grupos: \displaystyle SS_B = \sum_i \sum_j (\bar x_i - \bar x)^2

La parte clave del desarrollo matemático está en notar que la varianza total puede escribirse como la suma de las varianzas anteriores:
\displaystyle SS = \sum_i \sum_j (x_{i,j} - \bar x)^2 = SS_W + SS_B
En palabras, la agrupación de los datos explica parte de la varianza, por lo que una agrupación informativa explica una gran parte de la varianza. Nota: piensen en el análisis de varianza como una generalización al análisis de prueba t.
Definimos finalmente la estadística F como:
\displaystyle F = \frac{SS_B/(k-1)}{SS_W/(N-k)}
para N observaciones separadas en k grupos.

El comando en R que empleamos para el análisis de varianza es oneway.test. Veamos un ejemplo en donde tenemos observaciones separadas en 4 grupos (datos1, datos2, datos3, datos4). En este caso las observaciones son datos simulados empleando una distribución normal con $\latex \mu = 75, \sigma = 20$:

> datos1 = rnorm(20, mean=75, sd=20)
> datos2 = rnorm(20, mean=75, sd=20)
> datos3 = rnorm(20, mean=75, sd=20)
> datos4 = rnorm(20, mean=75, sd=20)
> 
> datos.tabla = data.frame(datos1, datos2, datos3, datos4)
> boxplot(datos.tabla)
> 
> datos.stack = stack(datos.tabla)
> names(datos.stack)
[1] "values" "ind"   
> 
> oneway.test(values ~ ind, data=datos.stack, var.equal=T)

        One-way analysis of means

data:  values and ind 
F = 0.9687, num df = 3, denom df = 76, p-value = 0.412

anova

Dado que p > 0.05 aceptamos la hipótesis nula que las medias de todos los grupos son iguales; en símbolos \bar x_1 = \bar x_2 = \cdots = \bar x_k.

Nota: Otros comandos que proporcionan resultados equivalentes a oneway.test son:

  • summary(aov(values ~ ind, data=datos.stack))
  • anova(lm(values ~ ind, data=datos.stack))
Tagged with: , ,

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s

A %d blogueros les gusta esto: