jueves, 6 de abril de 2017

estUNA

¿Qué es estUNA?


estUNA es una librería en R elaborada por mí para ampliar algunas funcionalidades del lenguaje R y facilitar la elaboración de los trabajos de estadística y el estudio independiente y a distancia de las materias de probabilidades y estadística de la Universidad Nacional Abierta en general. La librería contiene funciones para agrupar datos, generar resúmenes de datos agrupados y no agrupados, varios tipos de gráficas (histogramas, tortas, boxplots), contrastes de bondad de ajuste y de independencia chi-cuadrado, análisis de regresión, entre otras. La librería se actualiza cada semestre para incluir la data del trabajo práctico.

Eventualmente será publicada en el repositorio CRAN como un paquete. Actualmente, el archivo imagen (que permite trabajar con la librería) está disponible para su descarga en https://raw.githubusercontent.com/unamatematicaseltigre/estUNA/master/estUNA.

Para comenzar a trabajar con estUNA


Si ha instalado R en su computadora y desea trabajar sin conexión internet, debe haber previamente descargado el archivo imagen de estUNA desde la dirección http indicada arriba. Seguidamente debe mover este archivo al directorio de trabajo para que pueda acceder a él fácilmente desde la consola R (en Windows, el directorio de trabajo es usualmente la carpeta "Mis documentos"). Al abrir el programa , escriba el siguiente comando en la cónsola:
load("estUNA")

Esto cargará la librería estUNA desde el archivo imagen en el directorio de trabajo, para eventualmente poder trabajar con ella. Si surge algún error al efecto de "no se encontró el archivo", asegúrese que el archivo imagen de estUNA esté en el directorio de trabajo. Para salir de dudas, puede consultar cuál es el directorio de trabajo en R escribiendo el siguiente comando desde la consola R:
getwd()

Si está trabajando desde un RWeb server, la primera línea de código que debe colocar en el script a ser ejecutado es la siguiente:
load(url("https://raw.githubusercontent.com/unamatematicaseltigre/estUNA/master/estUNA"))

Recuerde que R es un lenguaje sensible a mayúsculas y minúsculas. Si intenta escribir los comandos de arriba cambiando las minúsculas por mayúsculas (o vice-versa), no obtendrá el resultado deseado. Tras cargar la librería estUNA mediante la instrucción load, necesita ubicar la data con la que se va a trabajar. La librería estUNA contiene la data utilizada para los trabajos prácticos de varios semestres. Cada uno de estos datasets es un dataframe cuyo identificador comienza por la letra "d" seguida del lapso académico en cuestión. Hasta ahora, los datasets disponibles en la librería son:
  • d20092- semestre 2009-2 de la UNA elaborado por Nivel Central (Ver enunciados).
  • d20102- semestre 2010-2 de la UNA elaborado por Nivel Central (Ver enunciados).
  • d20111- semestre 2011-1 de la UNA elaborado por mi para el Centro Local de Anzoategui. (Ver enunciados).
  • d20111b- semestre 2011-1 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20112- semestre 2011-2 de la UNA elaborado por mi para el Centro Local de Anzoategui. (Ver enunciados).
  • d20112b- semestre 2011-2 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20121- semestre 2012-1 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20131- semestre 2013-1 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20132- semestre 2013-2 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20141- semestre 2014-1 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20151e1 , d20151e2- semestre 2015-1 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20152- semestre 2015-2 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20161- semestre 2016-1 de la UNA elaborado por Nivel Central. (Ver enunciados).
  • d20171- semestre 2017-1 de la UNA elaborado por Nivel Central. (Ver enunciados).

En esta página ilustraremos algunas funcionalidades de la librería estUNA con la data del semestre 2010-2. Para los trabajos de este semestre, se trataba de determinar cuales factores inciden sobre las ventas de un vendedor (variable X1), basándose en una muestra de 60 vendedores de una compañía, incluyendo otras variables como la antigüedad de un vendedor(X2), su evaluación (X9) y la zona (X10) entre otras. Para referirnos fácilmente a las variables que componen el dataset de ese semestre, ejecutamos la instrucción attach (Ver en la página "Introducción al R"):
attach(d20102)

Estadística descriptiva sobre datos no agrupados


Cuando uno se enfrenta a un conjunto de datos por primera vez, es conveniente aplicar las técnicas de la estadística descriptiva antes de cualquier otro análisis posterior para "explorar" tendencias o comportamientos de esos datos y así poder formular hipótesis que tengan cierto asidero empírico. Una de esas técnicas es calcular las medidas estadísticas descriptivas que indican la tendencia central, la dispersión y la forma de esos datos. La librería estUNA provee una instrucción que genera todas estas medidas estadísticas en un cuadro resumen.

El resumen de los datos no agrupados se genera mediante la instrucción resumen. Así por ejemplo, si queremos generar un resumen con todas las medidas estadísticas relevantes para la variable X1, ejecutamos la instrucción resumen(X1), la cual arrojaría el siguiente resultado en la consola:

--------------------------------------------------------------
Resumen de la Variable X1  para datos sin agrupar

 Medidas de Tendencia Central 
  Moda    :  2752.402 
  Mediana :  2913.94 
  Media   :  3160.664 

 Medidas de Dispersion 
  Desv. Mediana Absoluta :  1048.406 
  Rango Inter Cuartil    :  1483.115 
  Desviacion Estandard   :  1384.611 

 Medidas de Posicion                                                                             
  Minimo     :
  421.98
  Cuartil 1  :
  2316.24
  Mediana    :  2913.94
  Cuartil 3  :  3799.355
  Maximo     :  6519.45

 Valores Atipicos
  Moderados :
    #53      #5      #9
6094,32 6125,96 6519,45
  Extremos  :
Ninguno

 Medidas de Forma
  Coeficiente de asimetria  :  0.548642
  Curtosis                  :  0.008867668

--------------------------------------------------------------

En el cuadro anterior se puede apreciar que aparte de reportar las medidas de tendencia central, dispersión y forma se indican los valores atípicos (moderados y extremos), si los hubiere. Los números con el prefijo # son los índices o números de renglón en el data frame correspondientes a los valores atípicos indicados bajo cada uno de ellos. Podemos observar también las diferencias entre la moda, la mediana y la media. Puesto que moda < mediana < media, se advierte la presencia de un sesgo a la derecha (Ver el libro de Webster, pg. 62). Esto se corrobora con el coeficiente de asimetría positivo, pero se puede apreciar mejor mediante una gráfica.

La instrucción graficar de estUNA genera varios tipos de gráficas, por lo menos las más relevantes para los cursos de estadística de pregrado. Aplicada a la variable X1, que es una variable cuantitativa, esta instrucción genera por defecto un histograma de frecuencias:
graficar(X1)


El histograma de frecuencias indica la ubicación de las tres medidas de tendencia central mediante líneas verticales de colores. Tal como aparece en la leyenda del gráfico, la línea roja sólida se corresponde a la moda, la verde a la mediana y la azul a la media. La curva roja punteada (----) se corresponde a la estimación de lo que se conoce como la densidad de núcleo (kernel density). Se supone que toda variable continua tiene una distribución de probabilidad continua y la densidad de núcleo es la estimación de la función de densidad subyacente de la variable. Se puede apreciar que la moda es el valor de la variable X1 para el cual la densidad de núcleo se hace máxima (observe el cruce entre la línea vertical roja y la curva roja punteada). Esto es un detalle importante: para las variables continuas, la moda no es el valor con mayor frecuencia (como el que calcula Excel), sino el valor para el cual la función de densidad (o su estimación, que es la densidad de núcleo), se hace máxima.

Por defecto, graficar, agrupa primero la variable en intervalos de clase según la regla de Sturges, que generalmente da mejores resultados. Sin embargo, es posible definir uno mismo su propio esquema de partición en intervalos de clase, inclusive usando intervalos de amplitud no uniforme. Por ejemplo, podríamos graficar la variable X1 agrupándola en dos intervalos de amplitud deferentes: [0,2000) y [2000,10000]. Esto se hace indicando un vector numérico en el argumento "intervalos":

graficar(X1,intervalos=c(0,2000,10000))


La instrucción graficar detecta hasta cierto punto si la variable a graficar es cualitativa o cuantitativa y según el caso, genera el tipo de gráfica más acorde con la escala de la variable. Por ejemplo, si graficamos la variable X10, que por ser la zona del vendedor es una variable categórica (nominal), la instrucción generaría una gráfica de torta:

graficar(X10)


A pesar de que la instrucción graficar intenta por defecto generar el tipo de gráfica más adecuado a la variable, es posible indicar uno mismo el tipo de gráfica para la variable. La instrucción graficar intentará, en la medida de lo posible, generar la gráfica deseada (Nota: no es posible generar histogramas para las variables categóricas y si deseamos generar una gráfica de torta para una variable cuantitativa, graficar intenta primero agrupar esa variable para "categorizarla").

graficar(X10,tipo="histograma")
Error en .local(x, ...) : Tipo de grafica invalido: histograma Calls: graficar -> graficar -> .local Ejecución interrumpida

Al querer graficar la variable X10 (tipo categórica) como un histograma de frecuencias se produce un error y se detiene la ejecución del script. Veamos que pasa si graficamos la variable X1 como una gráfica de torta:

graficar(X1,tipo="torta")

Como podemos observar, la gráfica de torta para la variable X1 es la menos adecuada para esta variable, porque no permite visualizar ni la tendencia central, ni la forma de la distribución, ni la dispersión. Existe otra modalidad de gráficas soportada por la instrucción graficar - las gráficas de barra. Este tipo de gráficas son más adecuadas para las variables ordinales. Si queremos generar una gráfica de barras para una variable x, tendríamos que ejecutar la instrucción graficar(x,tipo="barra").

Otro tipo importante de graficas son las gráficas de caja, o boxplots. Estas sirven más que todo para visualizar la dispersión y la atipicidad de los datos. Por ejemplo, para generar una gráfica de caja para la variable X1, ejecutamos la siguiente instrucción:
caja(X1)

En el diagrama de cajas, se indican los datos atipicos mediante puntos que yacen por encima o por debajo de los bigotes. Al lado de cada valor atipico señalado como un punto en la gráfica, se da el número de renglón en el data frame al cual corresponde el respectivo valor, para así poder facilitar la identificación de estos casos excepcionales.

Trabajando con datos agrupados


En la estadística descriptiva o inferencial, existen situaciones en las cuales es necesario agrupar los datos. Por ejemplo, para trabajar con los contrastes de bondad de ajuste o las pruebas de independencia basados en la distribución chi-cuadrado (para los estudiantes de la 746 o la 738/748), es preciso categorizar los datos primero, lo cual se hace agrupándolos. Por ejemplo, podemos agrupar la variable X1 y luego generar un resumen para datos agrupados:
X1_agrupada <- agrupar(X1)
resumen(X1_agrupada)
--------------------------------------------------------------
Resumen de la Variable X1 para datos agrupados.
Tabla de frecuencias:
        clase medio frec.abs   frec.rel
1  [0,1000]     500        2 0,03333333
2 (1000,2000]  1500        7 0,11666667
3 (2000,3000]  2500       24 0,40000000
4 (3000,4000]  3500       14 0,23333333
5 (4000,5000]  4500        5 0,08333333
6 (5000,6000]  5500        4 0,06666667
7 (6000,7000]  6500        4 0,06666667

 Medidas de Tendencia Central
  Moda        : 2629.63
  Mediana     : 2875
  Media       : 3183.333

 Medidas de Dispersion
  Rango Inter-cuartil  : 1607.143
  Desviacion Estandard : 1431.979

 Medidas de Posición
   Cuartil 1 : 2250
   Cuartil 3 : 3857.143
--------------------------------------------------------------

Obsérvese que primero se crea una variable identificada como X1_agrupada a la cual se le asigna el resultado de la función agrupar y luego se genera un resumen de la nueva variable X1_agrupada. Al igual que cuando generamos histogramas de frecuencia con un esquema de agrupación en intervalos de clase definida por nosotros, podemos definir nuestro propio esquema de agrupación y usarlo en la función agrupar. En lo que sigue daremos un ejemplo agrupando la variable X1 en cuatro intervalos de clase correspondientes a los 4 cuartiles, lo cual por cierto no tiene mucho sentido porque sabemos de antemano que cada clase tendría 1/4 de la frecuencia total:
otra_X1_agrupada <- agrupar(X1,intervalos=c(0,2316.24,2913.94,3799.355,10000))
resumen(otra_X1_agrupada)
--------------------------------------------------------------
Resumen de la Variable X1 para datos agrupados.
Tabla de frecuencias:
            clase   medio frec.abs frec.rel
1   [0,2316.2]    1158,10       15     0,25
2 (2316.2,2913.9] 2615,05       15     0,25
3 (2913.9,3799.4] 3356,65       15     0,25
4 (3799.4,10000]  6899,70       15     0,25

 Medidas de Tendencia Central 
  Moda        : 1158.1 2615.05 3356.65 6899.7 
  Mediana     : 2913.9 
  Media       : 3507.375 

 Medidas de Dispersion 
  Rango Inter-cuartil  : 1483.2 
  Desviacion Estandard : 2130.047 

 Medidas de Posición 
   Cuartil 1 : 2316.2                                                                            
   Cuartil 3 : 3799.4                                                                            
--------------------------------------------------------------  

Estadistica descriptiva para datos multivariantes


Cuando tenemos una muestra en la cual cada elemento tiene varias variables asociadas (lo que se conoce como datos multivariantes, seguramente estaremos interesados en estudiar el tipo de relación que se produce entre esas variables. En el marco de la estadística descriptiva, hay varias técnicas para abordar el estudio de las relaciones entre variables. La librería estUNA implementa un interfaz sencillo para dos de ellas: los diagramas de dispersión y los diagramas de caja comparativos.

Los diagramas de caja comparativos permiten estudiar la diferencia entre las distribuciones de una variable para grupos de individuos clasificados según otra variable. Por ejemplo, en el contexto de la data del semestre 2010-2, podriamos estar interesados en comparar las ventas de los vendedores (variable X1), según la zona del vendedor (variable X10). Para esto, ejecutamos la siguiente instrucción:
caja(X1,X10)

La instrucción caja(X1,X10) es la misma instrucción caja utilizada en el ejemplo univariante de arriba para generar el diagrama de caja de la variable X1, pero con un argumento adicional indicando una segunda variable de tipo categórica, como lo es la variable X10. Aquí también aparecen los valores atípicos, es decir, aquellas ventas consideradas como anómalas en cada zona. Esto ameritaría un análisis más profundo. Las lineas oscuras centrales representan las medianas, las hendiduras (o notches) son un indicativo de cuan diferentes son estas medianas. Según algunos autores (Chambers), si las hendiduras coinciden a la misma altura, las medianas no son significativamente distintas. Si bien la mediana de ventas para la zona 5 parece mucho mayor que la de las otras zonas, por el criterio de la coincidencia entre las hendiduras, esta diferencia no sería estadísticamente significativa, pues al parecer hay mucha variabilidad en las ventas de la zona 5.

También podemos indagar sobre la relación entre la antiguedad en años y las ventas (X1) mediante un diagrama de cajas comparativo. Sin embargo, como el segundo argumento de la instrucción caja tiene que ser una variable categorica, tenemos primero que "categorizar" la variable X2, que era una de las cosas que se pedía en el enunciado del trabajo de la 745 (semestre 2010-2). Esto se realiza mediante la instrucción agrupar. En el código siguiente, primero se transforma la variable X2, que es antigüedad en meses, a otra variable, que es antigüedad en años:
X2.anos <- X2 / 12
antiguedad <- agrupar(X2.anos,c(0,5,15,Inf))
caja(X1,antiguedad,titulox="Ventas",tituloa="Antiguedad")

En la última instrucción (caja(X1,antiguedad,titulox="Ventas",tituloa="Antiguedad")), se ilustra cómo colocar, de manera opcional, rotulos distintos de variables para el título de la gráfica. Los títulos en el los argumentos de la función son opcionales: titulox sería el nombre de la primera variable y tituloa el nombre de la segunda variable. Se puede observar que hay algún tipo de relación creciente entre la antigüedad y las ventas de un vendedor. Al parecer, las ventas para los vendedores con menos de 5 años de antigüedad son significativamente más bajas que las de los vendedores con mayor antigüedad. Aquí también observamos casos excepcionales (atípicos). ¿Qué ocurrirá por ejemplo con el vendedor #47, que tiene las ventas más bajas de todos pero mucha antigüedad? ¿Será que esta cansado y debe jubilarse?

Otra técnica útil para intuir relaciones entre variables es la gráfica de dispersión, aunque, a diferencia del diagrama de cajas comparativo, la gráfica de dispersión solo se aplica a variables continuas. La gráfica de dispersión se genera ubicando puntos en el plano cartesiano (X,Y), en donde los puntos (x,y) se corresponden a los valores observados de una variable X y otra variable Y para cada individuo de la muestra. Al ver la nube de puntos que resulta, se facilita la tarea de intuir una posible relación funcional entre ambas variables, lo cual sería de mucha utilidad para plantear modelos de regresión lineal. En las gráficas de dispersión generadas mediante la librería estUNA se incluye una curva de tendencia, la cual es una estimación de la tendencia central de la nube de puntos. Veamos un ejemplo de una gráfica de dispersión para las variables X1 (ventas) y X2 (antigüedad):
graficar.dispersion(X2,X1)
Con respecto a las técnicas gráficas expuestas aqui, es preciso acotar que su utilidad es principalmente "exploratoria". Los diagramas de caja comparativos y las gráficas de dispersión solo sirven para explorar o visualizar posibles relaciones entre variables. No demuestran si las relaciones son o no estadísticamente significativas. Para ello debemos recurrir a la regresión lineal o al análisis de varianza, que son técnicas de la inferencia estadística.

Regresión lineal con estUNA


La regresión lineal es una técnica que permite estimar los coeficientes de un modelo lineal propuesto. Un modelo lineal es una ecuación lineal que caracteriza la relación entre una variable (denominada variable dependiente) y otras variables, las variables dependientes. En una regresión lineal también se plantea la existencia de un error que explica de alguna manera las diferencias entre los valores observados de la variable dependiente y aquellos proyectados según el modelo lineal.
\[Y=\beta_0+\beta_1 X_1+\cdots+\beta_k X_k+\varepsilon\]
Además de estimar los coeficientes de un modelo lineal, la regresión lineal provee mecanismos para determinar si el modelo es adecuado, cuales variables son significativas y cuales no y si un modelo es mejor que otro. Vamos a ilustrar seguidamente el uso de estUNA para realizar una regresión lineal entre la variable X1 como variable dependiente y la variable X2 (antigüedad) como variable independiente:
modelo1 <- regresion.lineal(X1~X2)
resumen(modelo1)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo1 
  Marco de datos : variables globales 
  Formula        : X1 ~ X2 

Estimacion de los coeficientes poblacionales

              Estimacion Error Est. Estadistico T    p-valor
[Intercepto] 2419.231331 269.428216      8.979131 1.4355e-12
X2              6.647125   1.916005      3.469263 0.00099101

Prueba F global

  Valor F : 12.03579   gl. num: 1   gl. den : 58 
  p-valor : 0.00099101 

Coeficientes de determinacion
                                                                                                 
   R^2 : 0.171852   R^2 ajustado : 0.1575736                                                     
                                                                                                 
Residuos

  Minimo  : -3404.514 
  Mediana : -82.81871 
  Maximo  : 3434.928 

  Desv. estandar residual:  1270.849 

--------------------------------------------------------------

La expresión X1~X2 en el argumento de la función regresion.lineal se conoce como la fórumla del modelo lineal. En dicha fórmula, la variable a la izquierda del signo ~ es la variable dependiente, en este caso X1. La expresión a la derecha del signo ~ es una función lineal de las variables independientes. Por defecto, los modelos resultantes de la función regresion.lineal incluyen un coeficiente de intercepto. Para no incluir el coeficiente de intercepto, se debe agregar un -1 a la expresión derecha de la fórmula, la cual sería entonces X1~-1+X2. En este caso, el modelo1 es
\[X_1=\beta_0+\beta_1+X_2+\varepsilon\]
En este cuadro resumen se suministra la información más relevante con respecto a los resultados de la regresión lineal. En la columna "Estimacion" se indican los coeficientes estimados de las respectivas variables independientes. Según esta información, el modelo lineal estimado es:
\[X_1=2419.231331+6.647125 X_2\]
En las columnas a la derecha de la columna "Estimación" en esa misma tabla se suministra la información que permite hacer inferencia sobre los coeficientes poblacionales, a fin de determinar cuales variables son significativas y cuales no. De entre esta información, la más relevante es la referente al "p-valor". El p-valor indica cuán extremo es el valor estimado del estadístico T-Student en los contrastes de significatividad de los coeficientes poblacionales. Un p-valor muy bajo (por lo menos por debajo del nivel de significancia \(\alpha\) fijado por el investigador), indica que la estimación respectiva del estadístico T-Student es demasiado extrema, y por lo tanto debemos concluir que la variable independiente respectiva es significativa en el modelo. En este caso, ambos p-valores - \(1.4355\cdot 10^{-12}\) y \(0.00099101\) - indican que tanto el intercepto como la variable independiente X2 son signifivos. El p-valor de la prueba global F es muy bajo tambien, con lo cual reconfirmamos que las variables independientes de este modelo son significativas.
Podemos comparar estos resultados con la gráfica de dispersión de la variable X1 respecto a X2 de arriba. Observando los puntos de la nube de dispersión y la curva de tendencia verde en la región a la izquierda donde se encuentran los valores cercanos a cero de X2, constatamos que en efecto, para estos valores X1 es ligeramente mayor que 2000. Este es el intercepto de 2419.231331. Por otro lado, a medida que X2 aumenta, X1 aumenta ligeramente también. De hecho, por cada més de antiguedad adicional, el salario X1 aumentaría en aproximadamente 6 Bolívares. También podemos observar que la curva verde es aproximadamente lineal. Sin embargo, también podemos constatar que la dispersión de los puntos de la nube en torno a la "tendencia central" (la curva verde) es muy grande. ¿Cuanta variabilidad de X2 es explicada por la variabilidad de X1?

La respuesta a esta pregunta nos la da el coeficiente de determinación \(r^2\). Según los resultados en el cuadro resumen, aproximadamente un 17% de la variabilidad de X1 es explicada por la variabilidad de X2, lo cual indica que el modelo lineal estimado es de baja calidad (como regla general, un modelo lineal es bueno cuando el el coeficiente de determinación es superior al 70%). El coeficiente de determinación ajustado siempre es inferior al coeficiente de determinación. La diferencia entre ambos es mayor a medida que el modelo contiene más variables independientes, de modo que el coeficiente de determinación ajustado indica el grado de penalización impuesto por incluir muchas variables en un modelo. Un buen modelo de regresión es aquel para el cual el procentaje de variabilidad de la variable independiente explicado por las demás variables es alto, pero a la vez debe ser sencillo y fácil de interpretar. Los modelos que incluyen demasiadas variables independientes se tornan más complicados y por lo tanto su interpretación se hace más difícil.

Sin embargo, el modelo1 tiene muy pocas variables independientes y un coeficiente de determinación muy bajo, lo cual nos lleva a pensar que quizás debemos plantear un modelo lineal con más variables. Esto lo haremos seguidamente, pero primero realizaremos una tarea pendiente respecto a este modelo: el análisis de residuos.

La validez de los resultados en el cuadro resumen - los p-valores de cada coeficiente poblacional, el p-valor de la prueba global F y los coeficientes de determinación - dependen del grado en que los residuos del modelo se ajustan a ciertas suposiciones básicas que deben cumplir los modelos de regresión lineal. Por esta razón, ningún diagnóstico sobre un modelo de regresión lineal estaría completo si no tomamos en cuenta el análisis de residuos. En esta página no elaboraré sobre lo que conlleva el análisis de residuos. Para más información, invito al lector a consultar una pequeña monografía que he escrito sobre ello. En esta monografía se explica mediante varios ejemplos la importancia del análisis de residuos, las suposiciones que deben verificarse, las implicaciones de no cumplir con estas suposiciones y cómo realizar todo esto mediante la librería estUNA. Por los momentos, vamos a presentar a continuación las gráficas para el diagnóstico de los residuales del modelo1 generadas por la función graficar:

graficar(modelo1)
En las gráficas de diagnóstico de arriba, se puede apreciar cierta desviación de la normalidad de los residuos (por el p-valor bajo de 0,04). Sin embargo, podemos darnos cuenta de que esta desviación de la normalidad se debe a la presencia de residuales atípicos, los cuales ya habíamos observado en el diagrama de cajas comparativo de arriba. Por esta razón, no se considera la desviación de la normalidad de los residuos como particularmente severa y en todo caso, se pudiese mejorar la regresión extrayendo los individuos correspondientes a los residuos atipicos para estudiar estos aparte. En las graficas de residuos versus predicción y la variable X2, no se observa ninguna desviación importante de la homocedasticidad: la curva de tendencia verde es más o menos horizontal y coincide con el eje X de las gráficas y el ancho de la nube de puntos es uniforme a lo largo del eje X. Sin embargo, si se observan algunos residuales atipicos moderados más allá de las líneas punteadas amarillas, pero ya habiamos reparado sobre ellos en las gráficas anteriores. En conclusión, los residuos se ajustan a las suposiciones de un modelo de regresión lineal o por lo menos, las desviaciones no son importantes.

Como parte del análisis de residuos, deberiamos de revisar detenidamente aquellos individuos (vendedores) con residuos atípicos en el marco de la regesión lineal. Esto se hace mediante la función residuos.atipicos de estUNA, que devuleve un vector de índices del dataframe correspondiente a aquellos elementos con residuos atípicos. Este vector resultante se utiliza para indizar el data frame y así poder visualizar la data de estos vendedores:
d20102[residuos.atipicos(modelo1),]
        X1     X2       X3      X4    X5    X6     X7    X8   X9 X10 X11
5  6125.96 161.79 57805.10 7747.10  9.15  0.50 180.44 17.64 4.60   2   4
9  6519.45 127.64 46176.80 8846.20 12.54  1.24 203.25 17.42 4.90   4   3
34 5663.95 112.50 52866.64 6378.08 10.38  1.29 135.08 22.06 5.32   5   3
47  421.98 211.71 31649.65 6469.56  5.35 -0.22 191.30 14.56 3.56   4  17
53 6094.32  36.13 25179.04 2894.05  7.82  0.28  95.51 18.98 4.31   2  14

Se puede mejorar el modelo de regresión lineal en algo si incluimos otra variable- X7, que representa el número de cuentas asignadas al vendedor. Esto tiene sentido si consideramos que mientras más cuentas tenga un vendedor, mayor será su volúmen de ventas. Definimos el modelo 2:

modelo2 <- regresion.lineal(X1~1+X2+X7)
resumen(modelo2)
--------------------------------------------------------------
Resumen de regresion lineal

  MODELO         : modelo2 
  Marco de datos : variables globales 
  Formula        : X1 ~ 1 + X2 + X7 

Estimacion de los coeficientes poblacionales

              Estimacion Error Est. Estadistico T    p-valor
[Intercepto] 1631.195345 415.051536      3.930103 0.00023219
X2              4.603100   2.023349      2.274991 0.02668644
X7              8.475972   3.490704      2.428155 0.01835513

Prueba F global

  Valor F : 9.47385   gl. num: 2   gl. den : 57 
  p-valor : 0.00028039 

Coeficientes de determinacion

   R^2 : 0.2494835   R^2 ajustado : 0.2231496 

Residuos

  Minimo  : -3805.191 
  Mediana : -39.65316 
  Maximo  : 3487.275 

  Desv. estandar residual:  1220.384 

--------------------------------------------------------------

Se puede observar que, con un coeficiente de determinación algo superior, este modelo explica mejor la variabilidad de X1 que el primer modelo. Sin embargo, habría que realizar otro análisis de residuos.

Contrastes chi-cuadrado con estUNA: Bondad de ajuste


Las pruebas basadas en el estadístico chi-cuadrado de Pearson se usan, entre otras cosas, para:
  • Verificar si distribución de frecuencias de eventos respecto a una muestra aleatoria es consistente con cierta distribución teórica supuesta bajo la hipótesis nula. La hipótesis alternativa en este caso indicaría una discrepáncia significativa en la distribución de frecuencias muestrales respecto a la teórica. Estas son las denominadas pruebas de bondad de ajuste.
  • Realizar pruebas mediante una tabla de contingencia para verificar si las observaciones pareadas de dos variables categóricas son independientes entre sí.
Para el objetivo 8 del trabajo práctico del semestre 2010-2, se pedía verificar si la variable X1 sigue una distribución normal. El script a continuación ajusta X1 a una normal e imprime un resúmen y la información detallada sobre el ajuste (el objeto ajuste_X1):
ajuste_X1 <- ajustar(X1,"normal")
resumen(ajuste_X1)>br /> imprimir(ajuste_X1)
La salida del comando resumen(ajuste_X1)) muestra la información más relevante respecto al contraste de bondad de ajuste de la variable X1 a una distribución normal (otras distribuciones posibles son: "normal","expo" (la exponencial), "gamma", "unifc" (uniforme continua), "binom" (la binomial), "geom" (la geométrica) ,"pois" (la poisson) y "unifd" (uniforme discreta). En la tabla de resumen generada por este comando se indica la hipótesis nula (H_0), la hipótesis alternativa (H_1), el valor del estadístico chi-cuadrado, los grados de libertad del estadístico de pruebas y el p-valor.
--------------------------------------------------------------
Contraste de bondad de ajuste de X1 a una distribucion normal 
    
    H_0:  X1 ~ normal 
H_1:  X1 !~ normal 

Estdst. chi-cuadrado :  10.50914 
Grados de libertad   :  2 
p-valor              :  0.005223583 
--------------------------------------------------------------
Es muy instructivo imprimir un reporte detallado de los resultados del ajuste, para aprender sobre cómo funciona el test de bondad de ajuste chi-cuadrado. Para ajustar una muestra de alguna variable aleatoria a la distribución normal, es preciso primero estimar sus parámetros: la media y la desviación estándar. Estos se estiman a partir de la muestra de la manera usual (son estimadores de máximo-verosimilitud). Luego hay que dividir las observaciones en intervalos de clase - usualmente se utiliza la regla de Sturges para ello. Se cuentan las observaciones de la muestra en cada clase (estas son las frecuencias observadas) y se calcula la frecuencia esperada como el número de observaciones en cada clase que se "esperaría" tener según la hipótesis nula (en este caso, que la muestra proviene de una población normalmente distribuida). Sin embargo, hay un detalle- cuando alguna frecuencia esperada es menor que 5, se juntan los dos intervalos de clase contíguos en uno solo, pues las suposiciones básicas del contraste (respecto a que el estadístico de contraste es chi-cuadrado), se cumplen si todas las frecuencias esperadas de los intervalos de clases son mayores que cinco. Todo esto lo realiza la función ajustar() de estUNA. El reporte detallado para este ajuste se da a continuación:
--------------------------------------------------------------
Contraste de X1 a una distribucion normal 

Tabla de frecuencias observadas y esperadas:
        clase frec  frec.esp
1 (-Inf,2000]    9 12,056544
2 (2000,3000]   24 15,172192
3 (3000,4000]   14 16,439589
4 (4000,5000]    5 10,810430
5 (5000,Inf)     8  5,521245

Parametros:
   media     desv 
3160.664 1384.611 

Estadistico de prueba :  10.50914 
Grados de libertad    :  2 
p-valor               :  0.005223583 
--------------------------------------------------------------
Se dice que una imágen vale más de mil palabras y los contrastes de bondad de ajuste no son una excepción. Como cási cualquier objeto que provee la librería estUNA, podemos graficar estos contrastes:
graficar(ajuste_X1)
En la gráfica de arriba, la curva roja representa la función de densidad de una variable aleatoria normal con media y desviación estándar iguales a las de la muestra. Esta curva roja se coloca sobre el histograma generado a partir de la muestra de X1. Finalmente, en la parte superior se indican los parametros estimados y el p-valor.

Volvemos al asunto inicial: ¿es la distribución de X1 significativamente distinta a una normal? Todo depende del nivel de significancia con el cual estamos contrastando. Si nuestro nivel de significáncia es de 5%, cualquier p-valor menor a 5% sería considerado como evidencia significativa de una desviación de la distribución normal. Para el p-valor que obtuvimos (0,005223583), rechazariamos la hipótesis nula y concluiriamos que X1 no es normalmente distribuida. Las desviaciones de la normalidad se pueden apreciar en la gráfica de arriba: la distribución de X1 es sesgada y evidencia una asimetría hacia la derecha.

Contrastes chi-cuadrado con estUNA: Tests de Independencia


En la sección anterior mencionamos que los contrastes chi-cuadrado se utilizan también para determinar si dos variables categóricas son independientes, a través de las denominadas tablas de contingencia. Las filas de esta tabla representan los distintos niveles o categorías de una variable y las columnas las categorías de la otra variable. En cada celda anotamos la cuenta de la cantidad de observaciones que se corresponden a las respectivas categorías de la fila y la columna de esa celda. Bajo la hipótesis nula de este contraste, las variables categóricas son independientes y entonces esperariamos que cuenta esperada en cada celda fuese el producto de las frecuencias relativas marginales multiplicado por la cantidad total de datos (N). Como en el constraste de bondad de ajuste chi-cuadrado, se obtiene el estadístico de contraste sumando las diferencias cuadráticas de las frecuencias esperadas y observadas normalizadas por las respectivas frecuencias esperadas. Este proceso es tedioso de hacer a mano y por lo tanto la librería estUNA incluye una función para realizar estos contrastes: test.independencia.

En lo que sigue desarrollaremos en punto 8.2 del trabajo práctico de las asignaturas 738/748, en el cual se pedía "determinar si existe relación entre la Evaluación del vendedor (X9) y su Antigüedad (X2)". X9 y X2 son variables cuantitativas y por lo tanto, para poder aplicar el test de independencia chi-cuadrado, se deben transformar a variables categóricas. Seguidamente, en el enunciado se indicaba cómo debía realizarse la categorización:

Variable Categoría Rango
Evaluación (X9) Baja (1) x<3
Media (2) 3≤x<5
Alta (3) x≥5
Antigüedad (X2) Nivel de experiencia (1) Menor a 5 años
Nivel de experiencia (2) Entre 5 y 15 años
Nivel de experiencia (3) Más de 15 años

En el script a continuación, categorizamos las variables X2 y X9 de acuerdo al agrupamiento en la tabla y seguidamente desarrollamos el contraste de independiencia chi-cuadrado mediante la instrucción test.independencia. Como en los ejemplos anteriores, se supone que hemos cargado la librería estUNA (load("estUNA")) y que hemos indicado que trabajaremos con la data del semestre 2010-2 ( attach(d20102) ):

X2agr <- agrupar(X2,c(0,36,60,Inf))
X9agr <- agrupar(X9,c(0,3,5,Inf))
T1 <- test.independencia(X2agr,X9agr)
resumen(T1)
imprimir(T1)
El agrupamiento de la variable X2 se hace en meses, pues "meses" es la unidad en la que se expresa esta variable. Como la mayoría de los objetos de estUNA, los contrastes de independencia tienen su método de resúmen, el cual indica los resultados más importantes del contraste. Si queremos ver la tabla de frecuencias esperadas u otra información adicional, lo requerimos mediante la instrucción imprimir:
--------------------------------------------------------------
Contraste de independencia chi-cuadrado

H_0: X2 y X9 son independientes.
H_1: X2 y X9 no son independientes.
Estdst. chi-cuadrado : 4.072808 
Grados de libertad   : 4 
p-valor              : 0.3962421 
Observaciones:
(1) Hay frecuencias esperadas menores que 5.
--------------------------------------------------------------

--------------------------------------------------------------
Contraste de independencia chi-cuadrado

NOTA: En las tablas, las filas representan niveles de
    la variable X2 y las columnas representan los niveles
    de la variable  X9 .
     
Tabla de frecuencias observadas:
           
            [0,3] (3,5] (5, Inf)
  [0,36]     4     9     1      
  (36,60]    2     7     0      
  (60, Inf) 10    19     8      

Tabla de frecuencias esperadas:
          [0,3]     (3,5]     (5, Inf) 
[0,36]     3,733333  8,166667  2,100000
(36,60]    2,400000  5,250000  1,350000
(60, Inf)  9,866667 21,583333  5,550000

Estadistico de prueba :  4.072808 
Grados de libertad    :  4 
p-valor               :  0.3962421 
--------------------------------------------------------------

Según los resultados de este contraste, con un p-valor de 0,4 no hay evidencia significativa para rechazar la hipótesis nula. Por lo tanto, para efectos prácticos podemos afirmar que X2 y X9 son independientes. La instrucción test.independencia puede tomar como argumentos no solo variables agrupadas, sino también variables de tipo "factor", que en R se corresponden a las variables categóricas. Por ejemplo, si queremos constrastar la independencia entre la región geográfica del vendedor (variable X10) y la división dentro de cada región a la cual se corresponde el vendedor (X11), el contraste sería el siguiente:
T2 <- test.independencia(X10,X11)
resumen(T2)
--------------------------------------------------------------
Contraste de independencia chi-cuadrado

H_0: X10 y X11 son independientes.
H_1: X10 y X11 no son independientes.
Estdst. chi-cuadrado : 27.56548 
Grados de libertad   : 76 
p-valor              : 0.9999999 
Observaciones:
(1) Hay frecuencias esperadas menores que 5.
--------------------------------------------------------------

El ejemplo anterior es un poco tonto, pues determinar la independencia entre X10 y X11 no sería de mucho valor práctico. Sin embargo, sirve para ilustrar el uso del comando para variables de tipo "factor". Es de notar que el estadístico de contraste tiene muchos grados de libertad (76). Esto se debe a que X10 es una variable categorica con 5 niveles (las 5 zonas geográficas) y que X11 tiene 20 niveles (hasta 20 subdivisiones territoriales dentro de cada región). Por lo tanto, el estadístico de contraste tiene (5-1)x(20-1)=76 grados de libertad. Otro detalle es que el interprete R indica cuando alguna de las frecuencias esperadas es menor a 5, lo cual conlleva problemas como se explicó en la sección anterior sobre los contrastes de bondad de ajuste chi-cuadrado.

Envio de salida gráfica a archivos


A veces puede resultar conveniente enviar todos los gráficos a archivos. En R, esto se puede hacer abriendo un dispositivo gráfico, con lo cual se enviará cada página (o frame) de gráficos a archivos numerados (en el caso de archivos de extension .png o .jpg) o como una página dentro de un solo archivo en pdf. Los comandos de R para abrir dispositivos gráficos son pdf(...), png y jpeg(...)- los puntos suspensivos indican el nombre del archivo (o de los archivos) que desea crear. Tambien existen versiones de estos comandos para crear archivos gráficos de tipo bmp (windows bitmap) o tiff, pero estos dos generalmente no se usan.

Si queremos enviar la salida gráfica de nuestro script a archivos, debemos de ejecutar uno de los comandos mencionados arriba al principio de nuestro script. Al final, debemos cerrar el dispositivo gráfico (lo cual efectivamente escribe los archivos al disco duro) mediante la instrucción graphics.off(). Así por ejemplo, si queremos enviar todas las gráficas a un solo archivo pdf, esto se haría del siguiente modo (los puntos suspensivos entre ambas instrucciones indican el resto de nuestro script):

pdf("mi_archivo.pdf")
    :
    :
    :
graphics.off()
En el código de arriba, debe sustituir mi_archivo.pdf por un nombre de archivo válido de su escogencia, indicando la extensión de archivo pdf. Al terminar de ejecutar el script - y suponiendo que no contiene errores - el archivo pdf se creará en su directorio de trabajo. Los otros dispositivos gráficos no admiten multiples páginas como los archivos pdf. Por lo tanto, es necesario indicar que queremos una numeración consecutiva, lo cual se hace colocándo %02d o %03d (para numeraciones de 2 o 3 digitos respectivamente), como se muestra en el ejemplo siguiente, en el cual se envian todas las gráficas a archivos png:

png("trabajo%02d.png")
    :
    :
    :
graphics.off()
Esto facilita la elaboración del trabajo práctico, ya que desde el procesador de palabras que use, podrá insertar las gráficas desde archivos.

No hay comentarios:

Publicar un comentario