Introduccion A R y Python
Introduccion A R y Python
Índice General
Prólogo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1 Introducción y Preliminares . . . . . . . . . . . . . . . . 2
1.1 El entorno R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Programas relacionados. Documentación . . . . . . . . . . . . . . . . . . 2
1.3 Estadı́stica con R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.4 R en un sistema de ventanas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.5 Utilización interactiva de R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.6 Una sesión inicial. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.7 Ayuda sobre funciones y capacidades . . . . . . . . . . . . . . . . . . . . . 4
1.8 Órdenes de R. Mayúsculas y minúsculas. . . . . . . . . . . . . . . . . . . 5
1.9 Recuperación y corrección de órdenes previas . . . . . . . . . . . . . . 5
1.10 Ejecución de órdenes desde un archivo y redirección de la
salida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.11 Almacenamiento y eliminación de objetos . . . . . . . . . . . . . . . . 6
8 Distribuciones probabilı́sticas . . . . . . . . . . . . . . 37
8.1 Tablas estadı́sticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
8.2 Estudio de la distribución de unos datos . . . . . . . . . . . . . . . . . 38
8.3 Contrastes de una y de dos muestras . . . . . . . . . . . . . . . . . . . . . 41
iii
11 Modelos estadı́sticos en R . . . . . . . . . . . . . . . . 55
11.1 Definición de modelos estadı́sticos. Fórmulas . . . . . . . . . . . . 55
11.1.1 Contrastes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
11.2 Modelos lineales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
11.3 Funciones genéricas de extracción de información del modelo
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
11.4 Análisis de varianza. Comparación de modelos . . . . . . . . . . 60
11.4.1 Tablas ANOVA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
11.5 Actualización de modelos ajustados . . . . . . . . . . . . . . . . . . . . . 61
11.6 Modelos lineales generalizados . . . . . . . . . . . . . . . . . . . . . . . . . . 61
11.6.1 Familias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
11.6.2 La función glm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
11.7 Modelos de Mı́nimos cuadrados no lineales y de Máxima
verosimilitud . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
11.7.1 Mı́nimos cuadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
11.7.2 Máxima verosimilitud . . . . . . . . . . . . . . . . . . . . . . . . . 67
11.8 Algunos modelos no-estándar . . . . . . . . . . . . . . . . . . . . . . . . . . 67
iv
12 Procedimientos gráficos . . . . . . . . . . . . . . . . . . 69
12.1 Funciones gráficas de nivel alto . . . . . . . . . . . . . . . . . . . . . . . . . 69
12.1.1 La función plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
12.1.2 Representación de datos multivariantes . . . . . . . . . 70
12.1.3 Otras representaciones gráficas . . . . . . . . . . . . . . . . 70
12.1.4 Argumentos de las funciones gráficas de nivel alto
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
12.2 Funciones gráficas de nivel bajo . . . . . . . . . . . . . . . . . . . . . . . . 72
12.2.1 Anotaciones matemáticas . . . . . . . . . . . . . . . . . . . . . 74
12.2.2 Fuentes vectoriales Hershey. . . . . . . . . . . . . . . . . . . . 74
12.3 Funciones gráficas interactivas . . . . . . . . . . . . . . . . . . . . . . . . . . 74
12.4 Uso de parámetros gráficos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
12.4.1 Cambios permanentes. La función par() . . . . . . . 76
12.4.2 Cambios temporales. Argumentos de las funciones
gráficas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
12.5 Parámetros gráficos habituales . . . . . . . . . . . . . . . . . . . . . . . . . 76
12.5.1 Elementos gráficos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
12.5.2 Ejes y marcas de división. . . . . . . . . . . . . . . . . . . . . . 78
12.5.3 Márgenes de las figuras . . . . . . . . . . . . . . . . . . . . . . . 78
12.5.4 Figuras múltiples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
12.6 Dispositivos gráficos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
12.6.1 Inclusión de gráficos PostScript en documentos . . 81
12.6.2 Dispositivos gráficos múltiples . . . . . . . . . . . . . . . . . 82
12.7 Gráficos dinámicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
Apendice B Ejecución de R . . . . . . . . . . . . . . . . . 88
B.1 Ejecución de R en UNIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
B.2 Ejecución de R en Microsoft Windows . . . . . . . . . . . . . . . . . . . 91
Prólogo
Estas notas sobre R están escritas a partir de un conjunto de notas que describı́an los
entornos S y S-Plus escritas por Bill Venables y Dave Smith. Hemos realizado un pequeño
número de cambios para reflejar las diferencias entre R y S.
R es un proyecto vivo y sus capacidades no coinciden totalmente con las de S. En
estas notas hemos adoptado la convención de que cualquier caracterı́stica que se vaya a
implementar se especifica como tal en el comienzo de la sección en que la caracterı́stica es
descrita. Los usuarios pueden contribuir al proyecto implementando cualquiera de ellas.
Deseamos dar las gracias más efusivas a Bill Venables por permitir la distribución de
esta versión modificada de las notas y por ser un defensor de R desde su inicio.
Cualquier comentario o corrección serán siempre bienvenidos. Dirija cualquier correspon-
dencia a R-core@[Link].
Sugerencias al lector
La primera relación con R deberı́a ser la sesión inicial del Apendice A [Ejemplo de sesion],
página 84. Está escrita para que se pueda conseguir cierta familiaridad con el estilo de las
sesiones de R y para comprobar que coincide con la versión actual.
Muchos usuarios eligen R fundamentalmente por sus capacidades gráficas. Si ese es su
caso, deberı́a leer antes o después el Capı́tulo 12 [Graficos], página 69, sobre capacidades
gráficas y para ello no es necesario esperar a haber asimilado totalmente las secciones
precedentes.
Capı́tulo 1: Introducción y Preliminares 2
1 Introducción y Preliminares
1.1 El entorno R
R es un conjunto integrado de programas para manipulación de datos, cálculo y gráficos.
Entre otras caracterı́sticas dispone de:
• almacenamiento y manipulación efectiva de datos,
• operadores para cálculo sobre variables indexadas (Arrays), en particular matrices,
• una amplia, coherente e integrada colección de herramientas para análisis de datos,
• posibilidades gráficas para análisis de datos, que funcionan directamente sobre pantalla
o impresora, y
• un lenguaje de programación bien desarrollado, simple y efectivo, que incluye
condicionales, ciclos, funciones recursivas y posibilidad de entradas y salidas. (Debe
destacarse que muchas de las funciones suministradas con el sistema están escritas en
el lenguaje R)
El término “entorno” lo caracteriza como un sistema completamente diseñado y co-
herente, antes que como una agregación incremental de herramientas muy especı́ficas e
inflexibles, como ocurre frecuentemente con otros programas de análisis de datos.
R es en gran parte un vehı́culo para el desarrollo de nuevos métodos de análisis interactivo
de datos. Como tal es muy dinámico y las diferentes versiones no siempre son totalmente
compatibles con las anteriores. Algunos usuarios prefieren los cambios debido a los nuevos
métodos y tecnologı́a que los acompañan, a otros sin embargo les molesta ya que algún código
anterior deja de funcionar. Aunque R puede entenderse como un lenguaje de programación,
los programas escritos en R deben considerarse esencialmente efı́meros.
bibliotecas estándar) pero otras muchas están disponibles a través de Internet en CRAN
([Link]
Como hemos indicado, muchas técnicas estadı́sticas, desde las clásicas hasta la última
metodologı́a, están disponibles en R, pero los usuarios necesitarán estar dispuestos a traba-
jar un poco para poder encontrarlas.
Existe una diferencia fundamental en la filosofı́a que subyace en R (o S) y la de otros
sistemas estadı́sticos. En R, un análisis estadı́stico se realiza en una serie de pasos, con
unos resultados intermedios que se van almacenando en objetos, para ser observados o
analizados posteriormente, produciendo unas salidas mı́nimas. Sin embargo en SAS o SPSS
se obtendrı́a de modo inmediato una salida copiosa para cualquier análisis, por ejemplo,
una regresión o un análisis discriminante.
R preguntará si desea salvar los datos de esta sesión de trabajo. Puede responder yes
(Si), no (No) o cancel (cancelar) pulsando respectivamente las letras y, n o c, en cada
uno de cuyos casos, respectivamente, salvará los datos antes de terminar, terminará
sin salvar, o volverá a la sesión de R. Los datos que se salvan estarán disponibles en la
siguiente sesión de R.
Volver a trabajar con R es sencillo:
1. Haga que ‘trabajo’ sea su directorio de trabajo e inicie el programa como antes:
$ cd trabajo
$ R
2. Dé las órdenes que estime convenientes a R y termine la sesión con la orden q().
Bajo Microsoft Windows el procedimiento a seguir es básicamente el mismo: Cree una
carpeta o directorio. Ejecute R haciendo doble click en el icono correspondiente. Seleccione
New dentro del menú File para indicar que desea iniciar un nuevo problema (lo que eliminará
todos los objetos definidos dentro del espacio de trabajo) y a continuación seleccione Save
dentro del menú File para salvar esta imagen en el directorio que acaba de crear. Puede
comenzar ahora los análisis y cuando salga de R, éste le preguntará si desea salvar la imagen
en el directorio de trabajo.
Para continuar con este análisis posteriormente basta con pulsar en el icono de la imagen
salvada, o bien puede ejecutar R y utilizar la opción Open dentro del menú File para
seleccionar y abrir la imagen salvada.
> sink("[Link]")
enviará el resto de la salida, en vez de a la pantalla, al archivo del sistema operativo,
[Link], dentro del directorio de trabajo. La orden
> sink()
devuelve la salida de nuevo a la pantalla.
Si utiliza nombres absolutos de archivo en vez de nombres relativos, los resultados se
almacenaán en ellos, independientemente del directorio de trabajo.
2
El punto inicial del nombre de este archivo indica que es invisible en UNIX.
Capı́tulo 2: Cálculos sencillos. Números y vectores 7
1
Con argumentos diferentes, por ejemplo listas, la acción de c() puede ser diferente. Vea Sección 6.2.1
[Concatenacion de listas], página 29.
2
El sı́mbolo de subrayado, ‘_’, es un sinónimo del operador de asignación, pero no aconsejamos su uti-
lización ya que produce un código menos legible.
3
Aunque, de hecho, se almacena en .[Link] hasta que se ejecute otra orden.
Capı́tulo 2: Cálculos sencillos. Números y vectores 8
Advierta que max y min seleccionan el mayor y el menor valor de sus argumentos, incluso
aunque estos sean varios vectores. Las funciones paralelas pmax y pmin devuelven un vector
(de la misma longitud del argumento más largo) que contiene en cada elemento el mayor y
menor elemento de dicha posición de entre todos los vectores de entrada.
En la mayorı́a de los casos, el usuario no debe preocuparse de si los “números” de un vec-
tor numérico son enteros, reales o incluso complejos. Los cálculos se realizan internamente
como números de doble precisión, reales o complejos según el caso.
Para trabajar con números complejos, debe indicar explı́citamente la parte compleja.
Ası́
sqrt(-17)
devuelve el resultado NaN y un mensaje de advertencia, pero
sqrt(-17+0i)
realiza correctamente el cálculo de la raı́z cuadrada de este número complejo.
Capı́tulo 2: Cálculos sencillos. Números y vectores 9
Los operadores lógicos son < (menor), <= (menor o igual), > (mayor), >= (mayor o igual),
== (igual), y != (distinto). Además, si c1 y c2 son expresiones lógicas, entonces c1&c2 es
su intersección (“conjunción”), c1|c2 es su unión (“disyunción”) y !c1 es la negación de
c1.
Los vectores lógicos pueden utilizarse en expresiones aritméticas, en cuyo caso se trans-
forman primero en vectores numéricos, de tal modo que F se transforma en 0 y T en 1. Sin
embargo hay casos en que un vector lógico y su correspondiente numérico no son equiva-
lentes, como puede ver a continuación.
Por otra parte, la función paste() une todos los vectores de caracteres que se le sumi-
nistran y construye una sola cadena de caracteres. También admite argumentos numéricos,
que convierte inmediatamente en cadenas de caracteres. En su forma predeterminada, en
la cadena final, cada argumento original se separa del siguiente por un espacio en blanco,
aunque ello puede cambiarse utilizando el argumento sep="cadena", que sustituye el espacio
en blanco por cadena, la cual podrı́a ser incluso vacı́a.
Por ejemplo,
> labs <- paste(c("X","Y"), 1:10, sep="")
almacena, en labs, el vector de caracteres
c("X1", "Y2", "X3", "Y4", "X5", "Y6", "X7", "Y8", "X9", "Y10")
Recuerde que al tener c("X", "Y") solo dos elementos, deberá repetirse 5 veces para
obtener la longitud del vector 1:10.7
7
paste(..., collapse=ss) permite colapsar los argumentos en una sola cadena de caracteres separándolos
mediante ss. Además existen otras órdenes de manipulación de caracteres. como sub y substring. Puede
encontrar su descripción en la ayuda del programa.
Capı́tulo 2: Cálculos sencillos. Números y vectores 12
• Las listas son una forma generalizada de vector en las cuales los elementos no tienen
por qué ser del mismo tipo y a menudo son a su vez vectores o listas. Las listas
permiten devolver los resultados de los cálculos estadı́sticos de un modo conveniente.
Véase Sección 6.1 [Listas], página 28.
• Las hojas de datos (data frames) son estructuras similares a una matriz, en que cada
columna puede ser de un tipo distinto a las otras. Las hojas de datos son apropiadas
para describir ‘matrices de datos’ donde cada fila representa a un individuo y cada
columna una variable, cuyas variables pueden ser numéricas o categóricas. Muchos
experimentos se describen muy apropiadamente con hojas de datos: los tratamientos
son categóricos pero la respuesta es numérica. Véase Sección 6.3 [Hojas de datos],
página 29.
• Las funciones son también objetos de R que pueden almacenarse en el espacio de
trabajo, lo que permite extender las capacidades de R fácilmente. Véase Capı́tulo 10
[Escritura de funciones], página 46.
Capı́tulo 3: Objetos: Modos y atributos 14
modo, como para asignar un atributo a un objeto que carece de él. Es aconsejable consultar
la ayuda para familiarizarse con estas funciones.
y otras funciones genéricas, como summary(), producirán un resultado especial; todo ello
en función de la pertenencia a dicha clase.
Para eliminar temporalmente los efectos de la clase puede utilizar la función unclass().
Por ejemplo, si invierno pertenece a la clase "[Link]", entonces
> invierno
escribe el objeto en la forma de la clase, parecida a una matriz, en tanto que
> unclass(invierno)
lo imprime como una lista ordinaria. Sólo debe utilizar esta función en situaciones muy
concretas, como, por ejemplo, si hace pruebas para comprender el concepto de clase y de
función genérica.
Las clases y las funciones genéricas serán tratadas muy brevemente en la Sección 10.9
[Orientacion a objetos], página 54.
Capı́tulo 4: Factores Nominales y Ordinales 17
1
Para quienes no conocen la estructura administrativa de Australia, existen ocho estados y territorios
en la misma: Australian Capital Territory, New South Wales, Northern Territory, Queensland, South
Australia, Tasmania, Victoria, y Western Australia; y sus correspondientes abreviaturas son: act, nsw,
nt, qld, sa, tas, vic, y wa.
Capı́tulo 4: Factores Nominales y Ordinales 18
> MediaIngresos
act nsw nt qld sa tas vic wa
44.500 57.333 55.500 53.600 55.000 60.500 56.000 52.250
La función tapply() aplica una función, en este ejemplo la función mean(), a cada grupo
de componentes del primer argumento, en este ejemplo ingresos, definidos por los niveles
del segundo argumento, en este ejemplo FactorEstado, como si cada grupo fuese un vector
por sı́ solo. El resultado es una estructura cuya longitud es el número de niveles del factor.
Puede consultar la ayuda para obtener más detalles.
Suponga que ahora desea calcular las desviaciones tı́picas de las medias de ingresos por
estados. Para ello es necesario escribir una función en R que calcule la desviación tı́pica de
un vector. Aunque aún no se ha explicado en este texto cómo escribir funciones2 , puede
admitir que existe la función var() que calcula la varianza muestral o cuasi-varianza, y que
la función buscada puede construirse con la asignación:
> StdErr <- function(x) sqrt(var(x)/length(x))
Ahora puede calcular los valores buscados mediante
> ErrorTipicoIngresos <- tapply(ingresos, FactorEstado, StdErr)
con el siguiente resultado:
> ErrorTipicoIngresos
act nsw nt qld sa tas vic wa
1.500000 4.310195 4.500000 4.106093 2.738613 0.500000 5.244044 2.657536
Como ejercicio puede calcular el intervalo de confianza al 95% de la media de ingresos por
estados. Para ello puede utilizar la función tapply(), la función length() para calcular los
tamaños muestrales, y la función qt() para encontrar los percentiles de las distribuciones t
de Student correspondientes.
La función tapply() puede utilizarse para aplicar una función a un vector indexado por
diferentes categorı́as simultáneamente. Por ejemplo, para dividir la muestra tanto por el
estado como por el sexo. Los elementos del vector se dividirán en grupos correspondientes
a las distintas categorı́as y se aplicará la función a cada uno de dichos grupos. El resultado
es una variable indexada etiquetada con los niveles de cada categorı́a.
La combinación de un vector3 con un factor para etiquetarlo, es un ejemplo de lo que se
llama variable indexada desastrada (ragged array) puesto que los tamaños de las subclases
son posiblemente irregulares. Cuando estos tamaños son iguales la indexación puede hacerse
implı́citamente y además más eficientemente, como veremos a continuación.
creados por la función ordered() los denominaremos ordinales. En la mayorı́a de los casos
la única diferencia entre ambos tipos de factores consiste en que los ordinales se imprimen
indicando el orden de los niveles. Sin embargo los contrastes generados por los dos tipos de
factores al ajustar Modelos lineales, son diferentes.
Capı́tulo 5: Variables indexadas. Matrices 20
4
b es esencialmente el resultado del operador “barra hacia atrás” de Matlab.
Capı́tulo 5: Variables indexadas. Matrices 27
6.1 Listas
En R, una lista es un objeto consistente en una colección ordenada de objetos, conocidos
como componentes.
No es necesario que los componentes sean del mismo modo, ası́ una lista puede estar
compuesta de, por ejemplo, un vector numérico, un valor lógico, una matriz y una función.
El siguiente es un ejemplo de una lista:
> Lst <- list(nombre="Pedro", esposa="Marı́a", [Link]=3,
[Link]=c(4,7,9))
Los componentes siempre están numerados y pueden ser referidos por dicho número.
En este ejemplo, Lst es el nombre de una lista con cuatro componentes, cada uno de los
cuales puede ser referido, respectivamente, por Lst[[1]], Lst[[2]], Lst[[3]] y Lst[[4]].
Como, además, Lst[[4]] es un vector, Lst[[4]][1] refiere su primer elemento.
La función length() aplicada a una lista devuelve el número de componentes (del primer
nivel) de la lista.
Los componentes de una lista pueden tener nombre, en cuyo caso pueden ser referidos
también por dicho nombre, mediante una expresión de la forma
nombre de lista$nombre de componente
Esta convención permite la obtención de una componente sin tener que recurrir a su
número.
En el ejemplo anterior,
Lst$nombre coincide con Lst[[1]] y vale "Pedro",
Lst$esposa coincide con Lst[[2]] y vale "Marı́a",
Lst$[Link][1] coincide con Lst[[4]][1] y vale 4.
También es posible utilizar los nombres de los componentes entre dobles corchetes, por
ejemplo, Lst[["nombre"]] coincide con Lst$nombre. Esta opción es muy útil en el caso
en que el nombre de los componentes se almacena en otra variable, como en
> x <- "nombre"; Lst[[x]]
Es muy importante distinguir claramente entre Lst[[1]] y Lst[1]. ‘[[. . . ]]’ es el
operador utilizado para seleccionar un sólo elemento, mientras que ‘[. . . ]’ es un operador
general de indexado. Esto es, Lst[[1]] es el primer objeto de la lista Lst, y si es una lista
con nombres, el nombre no está incluido. Por su parte, Lst[1], es una sublista de la lista
Lst consistente en la primera componente. Si la lista tiene nombre, éste se transfiere a la
sublista.
Los nombres de los componentes pueden abreviarse hasta el mı́nimo de letras necesarios
para identificarlos de modo exacto. Ası́, en
> Lista <- list(coeficientes=c(1.3,4), covarianza=.87)
Lst$coeficientes puede especificarse mediante Lista$coe, y Lista$covarianza como
Lista$cov.
El vector de nombres es un atributo de la lista, y como el resto de atributos puede ser
manipulado. Además de las listas, también otras estructuras pueden poseer el atributo
names.
Capı́tulo 6: Listas y hojas de datos 29
1
Aunque R permite lo contrario, deberı́an ser distintos entre sı́
2
Hemos utilizado esta traducción por analogı́a con la “hoja de cálculo”
Capı́tulo 6: Listas y hojas de datos 30
y, como vimos, ls (o objects) puede usarse para examinar los contenidos de cualquier
posición en la trayectoria de búsqueda.
Por último, desconectamos la hoja de datos y comprobamos que ha sido eliminada de la
trayectoria de búsqueda.
> detach("lentejas")
> search()
[1] ".GlobalEnv" "Autoloads" "package:base"
Capı́tulo 7: Lectura de datos de un archivo 33
Para acceder a los datos incluidos en una biblioteca, basta utilizar el argumento package
en la función data. Por ejemplo,
data(package="nls")
data(Puromycin, package="nls")
Si la biblioteca ya ha sido conectada mediante la función library, sus datos habrán
sido incluidos automáticamente en la trayectoria de búsqueda y no será necesario incluir el
argumento package. Ası́,
library(nls)
data()
data(Puromycin)
presentará una lista de todos los datos de todas las bibliotecas conectadas en ese momento
(que serán al menos la biblioteca base y la biblioteca nls) y posteriormente cargará los datos
Puromycin de la primera librerı́a en la trayectoria de búsqueda en que encuentre unos datos
con dicho nombre.
Las librerı́as creadas por los usuarios son una fuente valiosa de datos. Por supuesto, las
notas del Dr. Venables, fuente original de esta introducción, contienen un conjunto de datos
que se encuentra disponible en CRAN en la biblioteca Rnotes.
2
En S-Plus la carga explı́cita no es necesaria.
Capı́tulo 7: Lectura de datos de un archivo 36
3
Acrónimo en inglés de Standard Query Language.
Capı́tulo 8: Distribuciones probabilı́sticas 37
8 Distribuciones probabilı́sticas
> qf(0.99, 2, 7)
16 | 070355555588
18 | 000022233333335577777777888822335777888
20 | 00002223378800035778
22 | 0002335578023578
24 | 00228
26 | 23
28 | 080
30 | 7
32 | 2337
34 | 250077
36 | 0000823577
38 | 2333335582225577
40 | 0000003357788888002233555577778
42 | 03335555778800233333555577778
44 | 02222335557780000000023333357778888
46 | 0000233357700000023578
48 | 00000022335800333
50 | 0370
En vez del diagrama de tallo y hojas, puede representar el histograma utilizando la
función hist.
> hist(eruptions)
# define intervalos menores y a~
nade un gráfico de densidad
> hist(eruptions, seq(1.6, 5.2, 0.2), prob=TRUE)
> lines(density(eruptions, bw=0.1))
> rug(eruptions) # muestra los puntos
La función density permite realizar gráficos de densidad y la hemos utilizado para
superponer este gráfico en el ejemplo. La anchura de banda, bw, ha sido elegida probando
varias, ya que el valor predeterminado produce un gráfico mucho más suavizado. Si necesita
Capı́tulo 8: Distribuciones probabilı́sticas 39
utilizar métodos automáticos de elección de ancho de banda, utilice las bibliotecas MASS
y KernSmooth.
Histogram of eruptions
0.7
0.6
0.5
Relative Frequency
0.4
0.3
0.2
0.1
0.0
eruptions
0.4
0.2
0.0
Los gráficos cuantil-cuantil ("Q-Q plots") pueden ayudarnos a examinar los datos más
cuidadosamente.
Capı́tulo 8: Distribuciones probabilı́sticas 40
> par(pty="s")
> qqnorm(long); qqline(long)
que muestran un ajuste razonable, aunque la cola de la derecha es más corta de lo que serı́a
esperable en una distribución normal. Vamos a compararla con unos datos pseudoaleatorios
tomados de una distribución t5 .
Normal Q−Q Plot
5.0
4.5
Sample Quantiles
4.0
3.5
3.0
−2 −1 0 1 2
Theoretical Quantiles
data: long
W = 0.9793, p-value = 0.01052
y el contraste de Kolmogorov-Smirnov
> [Link](long, "pnorm", mean=mean(long), sd=sqrt(var(long)))
data: long
D = 0.0661, p-value = 0.4284
alternative hypothesis: [Link]
Capı́tulo 8: Distribuciones probabilı́sticas 41
Hemos utilizado los datos como ejemplo de uso de las funciones, sin estudiar si el mismo es
válido. En este caso no lo serı́a ya que se han estimado los parámetros de la distribución
normal a partir de la misma muestra.
> boxplot(A, B)
que muestra claramente que el primer grupo tiende a tener mayores resultados que el se-
gundo.
80.04
80.02
80.00
79.98
79.96
79.94
1 2
Para contrastar la igualdad de medias de las dos poblaciones, se puede utilizar el con-
traste t de Student para dos muestras independientes, del siguiente modo:
> [Link](A, B)
data: A and B
t = 3.2499, df = 12.027, p-value = 0.00694
Capı́tulo 8: Distribuciones probabilı́sticas 42
data: A and B
F = 0.5837, num df = 12, denom df = 7, p-value = 0.3938
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.1251097 2.1052687
sample estimates:
ratio of variances
0.5837405
que no muestra evidencia de diferencias significativas (Bajo las condiciones del modelo, que
incluyen normalidad). Si hubiésemos admitido esta hipótesis previamente, podrı́amos haber
realizado un contraste más potente, como el siguiente:
> [Link](A, B, [Link]=TRUE)
data: A and B
t = 3.4722, df = 19, p-value = 0.002551
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.01669058 0.06734788
sample estimates:
mean of x mean of y
80.02077 79.97875
Como hemos indicado, una de las condiciones de aplicación de los contrastes anteriores
es la normalidad. Si ésta falla, puede utilizar el contraste de dos muestras de Wilcoxon (o
de Mann-Whitney) que solo presupone en la hipótesis nula que la distribución común es
continua.
> library(ctest)
> [Link](A, B)
data: A and B
W = 89, p-value = 0.007497
alternative hypothesis: true mu is not equal to 0
Warning message:
Cannot compute exact p-value with ties in: [Link](A, B)
Advierta el mensaje de advertencia (Warning . . . ): Existen valores repetidos en cada mues-
tra, lo que sugiere que los datos no proceden de una distribución continua (puede que ello
ocurra debido al redondeo).
Además del diagrama de cajas, existen más métodos para comparar gráficamente dos
muestras. Ası́, las órdenes siguientes:
> library(stepfun)
> plot(ecdf(A), [Link]=FALSE, verticals=TRUE, xlim=range(A, B))
> plot(ecdf(B), [Link]=FALSE, verticals=TRUE, add=TRUE)
representan las dos funciones de distribución empı́ricas. Por otra parte la función qqplot
realizarı́a un gráfico cuantil-cuantil de las dos muestras.
El contraste de Kolmogorov-Smirnov, que sólo presupone que la distribución común es
continua, también puede aplicarse:
> [Link](A, B)
data: A and B
D = 0.5962, p-value = 0.05919
alternative hypothesis: [Link]
Warning message:
cannot compute correct p-values with ties in: [Link](A, B)
siéndole de aplicación la misma precaución del contraste de Wilcoxon.
Capı́tulo 9: Ciclos. Ejecución condicional 44
2
No existe equivalente para esta orden en Fortran o Basic
Capı́tulo 10: Escritura de nuevas funciones 46
E = Iv − R−1/2 N 0 K −1 N R−1/2 = Iv − A0 A,
donde A = K −1/2 N R−1/2 .
Por ejemplo, la función podrı́a escribirse ası́:
> EfiDisBlo <- function(bloques, variedades) {
bloques <- [Link](bloques) # pequeña precaución
b <- length(levels(bloques))
variedades <- [Link](variedades) # pequeña precaución
v <- length(levels(variedades))
K <- [Link](table(bloques)) # elimina el atributo dim
R <- [Link](table(variedades)) # elimina el atributo dim
N <- table(bloques, variedades)
A <- 1/sqrt(K) * N * rep(1/sqrt(R), rep(b, v))
sv <- svd(A)
list(eficiencia=1 - sv$d^2, cvbloques=sv$u, cvvariedad=sv$v)
}
Desde el punto de vista numérico, es levemente mejor trabajar con la función descom-
posición SVD en vez de con la función de los autovalores.
El resultado de esta función es una lista que contiene los factores de eficiencia como
primera componente, y que además incluye dos contrastes, puesto que, a veces, suministran
información adicional útil.
tiempo que ilustra el hecho de que las funciones pueden ser cortas y al mismo tiempo muy
efectivas y útiles.
SinNombres <- function(a) {
## Elimina los nombres de dimensiones para impresión compacta.
d <- list()
l <- 0
for(i in dim(a)) {
d[[l <- l + 1]] <- rep("", i)
}
dimnames(a) <- d
a
}
Una vez definida la función, para imprimir la matriz X en forma compacta basta con
escribir
> SinNombres(X)
Esta función es de especial utilidad al imprimir variables indexadas de tipo entero y de
gran tamaño, en que el interés real se centra más en los posibles patrones que en los valores
en sı́ mismos.
10.7 Ámbito
Este apartado es algo más técnico que otras partes de este documento. Sin embargo,
pormenoriza una de las mayores diferencias entre S-Plus y R.
Los sı́mbolos que tienen lugar en el cuerpo de una función se dividen en tres clases:
parámetros formales, variables locales y variables libres. Los parámetros formales son los
que aparecen en la lista de argumentos de la función y sus valores quedan determinados
por el proceso de asignación de los argumentos de la función a los parámetros formales.
Las variables locales son aquellas cuyos valores están determinados por la evaluación de
expresiones en el cuerpo de las funciones. Las variables que no son parámetros formales
ni variables locales se denominan variables libres. Las variables libres se transforman en
variables locales si se les asigna valor. Para aclarar los conceptos, consideremos la siguiente
función:
f <- function(x) {
y <- 2*x
print(x)
print(y)
print(z)
}
En esta función, x es un parámetro formal, y es una variable local y z es una variable
libre.
En R la asignación de valor a una variable libre se realiza consultando el entorno en el
que la función se ha creado, lo que se denomina ámbito léxico. En primer lugar definamos
la función cubo:
cubo <- function(n) {
sq <- function() n*n
n*sq()
}
La variable n de la función sq no es un argumento para esta función. Por tanto es una
variable libre y las reglas de ámbito deben utilizarse para determinar el valor asociado con
ella. En un ámbito estático (como en S-Plus) el valor es el asociado con una variable global
llamada n. En un ámbito léxico (como en R) es un parámetro para la función cubo puesto
que hay una asignación activa para la variable n en el momento en que se define la función
sq. La diferencia de evaluación entre R y S-Plus es que S-Plus intenta encontrar una
variable global llamada n en tanto que R primero intenta encontrar una variable llamada n
en el entorno creado cuando se activó cubo.
## primera evaluación en S
Capı́tulo 10: Escritura de nuevas funciones 52
S> cubo(2)
Error in sq(): Object "n" not found
Dumped
S> n <- 3
S> cubo(2)
[1] 18
## la misma función evaluada en R
R> cubo(2)
[1] 8
El ámbito lexicográfico puede utilizarse para conceder a las funciones un estado cam-
biante. En el siguiente ejemplo mostramos cómo puede utilizarse R para simular una cuenta
bancaria. Una cuenta bancaria necesita tener un balance o total, una función para realizar
depósitos, otra para retirar fondos, y una última para conocer el balance.
Conseguiremos esta capacidad creando tres funciones dentro de [Link] y de-
volviendo una lista que los contiene. Cuando se ejecuta [Link] toma un argumento
numérico, total, y devuelve una lista que contiene las tres funciones. Puesto que estas fun-
ciones están definidas dentro de un entorno que contiene a total, éstas tendrán acceso a
su valor.
El operador de asignación especial, <<-, se utiliza para cambiar el valor asociado con
total. Este operador comprueba los entornos creados desde el actual hasta el primero
hasta encontrar uno que contenga el sı́mbolo total y cuando lo encuentra, sustituye su
valor en dicho entorno por el valor de la derecha de la expresión. Si se alcanza el nivel
superior, correspondiente al entorno global, sin encontrar dicho sı́mbolo, entonces lo crea
en él y realiza la asignación. Para muchos usos, <<- crea una variable global y le asigna
el valor de la derecha de la expresión4 . Solo cuando <<- ha sido utilizado en una función
que ha sido devuelta como el valor de otra función ocurrirá la conducta especial que hemos
descrito.
[Link] <- function(total) {
list(
deposito = function(importe) {
if(importe <= 0)
stop("Los depósitos deben ser positivos!\n")
total <<- total + importe
cat("Depositado",importe,". El total es", total, "\n\n")
},
retirada = function(importe) {
if(importe > total)
stop("No tiene tanto dinero!\n")
total <<- total - importe
cat("Descontado", importe,". El total es", total,"\n\n")
},
balance = function() {
cat("El total es", total,"\n\n")
}
)
4
En cierto sentido esto emula la conducta en S-Plus puesto que en S-Plus este operador siempre crea o
asigna a una variable global.
Capı́tulo 10: Escritura de nuevas funciones 53
Antonio$retirada(30)
Antonio$balance()
Roberto$balance()
Antonio$deposit(50)
Antonio$balance()
Antonio$retirada(500)
11 Modelos estadı́sticos en R
En este apartado, suponemos al lector familiarizado con la terminologı́a estadı́stica, en
particular con el análisis de regresión y el análisis de varianza. Posteriormente haremos
algunas suposiciones más ambiciosas, particularmente el conocimiento de modelos lineales
generalizados y regresión no lineal.
Los requisitos para el ajuste de modelos estadı́sticos están suficientemente bien definidos
para hacer posible la construcción de herramientas generales de aplicación a un amplio
espectro de problemas.
R contiene un conjunto de posibilidades que hace que el ajuste de modelos estadı́sticos
sea muy simple. Como hemos mencionado en la introducción, la salida básica es mı́nima, y
es necesario utilizar las funciones extractoras para obtener todos los detalles.
y = Xβ + e
Ejemplos
Antes de dar una definición formal, algunos ejemplos ayudarán a centrar las ideas.
Supongamos que y, x, x0, x1, x2, . . . son variables numéricas, que X es una matriz y que
A, B, C, . . . son factores. Las fórmulas que aparecen en la parte izquierda de la siguiente
tabla, especifican los modelos estadı́sticos descritos en la parte de la derecha.
y∼x
y∼1+x Ambos definen el mismo modelo de regresión lineal de y sobre x. El primero
contiene el término independiente implı́cito y el segundo, explı́cito.
y∼0+x
y ∼ -1 + x
y ∼ x - 1 Regresión lineal de y sobre x sin término independiente, esto es, que pasa por
el origen de coordenadas.
log(y) ∼ x1 + x2
Regresión múltiple de la variable transformada, log(y), sobre x 1 y x 2 (con un
término independiente implı́cito).
Capı́tulo 11: Modelos estadı́sticos en R 56
y ∼ poly(x,2)
y ∼ 1 + x + I(x^2)
Regresión polinomial de y sobre x de segundo grado. La primera forma utiliza
polinomios ortogonales y la segunda utiliza potencias de modo explı́cito.
y ∼ X + poly(x,2)
Regresión múltiple de y con un modelo matricial consistente en la matriz X y
términos polinomiales en x de segundo grado.
y∼A Análisis de varianza de entrada simple de y, con clases determinadas por A.
y∼A+x Análisis de covarianza de entrada simple de y, con clases determinadas por A,
y con covariante x.
y ∼ A*B
y ∼ A + B + A:B
y ∼ B %in% A
y ∼ A/B Modelo no aditivo de dos factores de y sobre A y B. Los dos primeros es-
pecifican la misma clasificación cruzada y los dos últimos especifican la misma
clasificación anidada. En términos abstractos, los cuatro especifican el mismo
subespacio de modelos.
y ∼ (A + B + C)^2
y ∼ A*B*C - A:B:C
Experimento con tres factores con un modelo que contiene efectos principales e
interacciones de dos factores solamente. Ambas fórmulas especifican el mismo
modelo.
y∼A*x
y ∼ A/x
y ∼ A/(1 + x) - 1
Modelos de regresión lineal simple separados de y sobre x para cada nivel de
A. La última forma produce estimaciones explı́citas de tantos términos inde-
pendientes y pendientes como niveles tiene A.
y ∼ A*B + Error(C)
Un experimento con dos factores de tratamiento, A y B, y estratos de error
determinados por el factor C. Por ejemplo, un experimento split plot, con
gráficos completos (y por tanto también subgráficos) determinados por el factor
C.
El operador ∼ se utiliza para definir una fórmula de modelo en R. La forma, para un
modelo lineal ordinario es
respuesta ∼ op 1 term 1 op 2 term 2 op 3 term 3 . . .
donde
respuesta es un vector o una matriz (o una expresión que evalúe a un vector o matriz)
que definen, respectivamente, la o las variables respuesta
op i es un operador, bien +, bien -, que implica la inclusión o exclusión, respectiva-
mente, de un término en el modelo. El primero, +, es opcional.
term i es un término de uno de los siguientes tipos
Capı́tulo 11: Modelos estadı́sticos en R 57
11.1.1 Contrastes
Es necesario conocer, aunque sea someramente, el modo en que las fórmulas del modelo
determinan las columnas de la matriz del modelo. Esto es sencillo si las variables son
Capı́tulo 11: Modelos estadı́sticos en R 58
continuas, ya que cada una constituirá una columna de dicha matriz. Del mismo modo, si
el modelo incluye un término independiente, contribuirá con una columna de unos.
En el caso de un factor, A, con k niveles, la respuesta depende de si el factor es nominal
u ordinal. En el caso de un factor nominal, se generan k − 1 columnas correspondien-
tes a los indicadores desde el segundo hasta el k-ésimo nivel del factor. (Por tanto, la
parametrización implı́cita consiste en contrastar la respuesta del primer nivel frente a cada
uno de los restantes niveles.) En el caso de un factor ordinal, las k − 1 columnas son los
polinomios ortogonales sobre 1, ..., k, omitiendo el término constante.
Esta situación puede parecerle complicada, pero aún hay más. En primer lugar, si el
término independiente se omite en un modelo que contiene algún término de tipo factor, el
primero de dichos términos se codifica en k columnas correspondientes a los indicadores de
todos los niveles del factor. En segundo lugar, todo este comportamiento puede cambiarse
mediante el argumento contrasts de options. Los valores predeterminados son:
options(contrasts = c("[Link]", "[Link]"))
La razón por la que se indican estos valores es que los valores predeterminados en R son
distintos de los de S en el caso de factores nominales, ya que S utiliza los contrastes de
Helmert. Por tanto, para obtener los mismos resultados que en S-Plus, deberá escribir:
options(contrasts = c("[Link]", "[Link]"))
Esta diferencia es deliberada, ya que entendemos que los contrastes predeterminados de R
son más sencillos de interpretar para los principiantes.
Caben aún más posibilidades, ya que el esquema de contraste a utilizar puede fijarse
para cada término del modelo utilizando las funciones contrasts y C.
Tampoco hemos considerado los términos de interacción, que generan los productos de
las columnas introducidas por los términos de sus componentes.
Pese a que los detalles son complicados, las fórmulas de modelos en R generan habi-
tualmente los modelos que un estadı́stico experto podrı́a esperar, supuesto que se preserve
la marginalidad. Por ejemplo, el ajuste de un modelo con interacción y, sin embargo, sin
los correspondientes efectos principales conducirá en general a resultados sorprendentes, y
debe reservarse sólo a los especialistas.
anova(objeto 1, objeto 2)
Compara un submodelo con un modelo externo y produce una tabla de análisis
de la varianza.
coefficients(objeto)
Extrae la matriz de coeficientes de regresión.
Forma reducida: coef(objeto).
deviance(objeto)
Suma de cuadrados residual, ponderada si es lo apropiado.
formula(objeto)
Extrae la fórmula del modelo.
plot(objeto)
Crea cuatro gráficos que muestran los residuos, los valores ajustados y algunos
diagnósticos
predict(objeto, newdata=[Link])
La nueva hoja de datos que se indica debe tener variables cuyas etiquetas coin-
cidan con las de la original. El resultado es un vector o matriz de valores
predichos correspondiente a los valores de las variables de [Link].
print(objeto)
Imprime una versión concisa del objeto. A menudo se utiliza implı́citamente.
residuals(objeto)
Extrae la matriz de residuos, ponderada si es necesario.
La forma reducida es resid(objeto).
step(objeto)
Selecciona un modelo apropiado añadiendo o eliminando términos y preservando
las jerarquı́as. Se devuelve el modelo que en este proceso tiene el máximo valor
de AIC1 .
summary(objeto)
Imprime un resumen estadı́stico completo de los resultados del análisis de re-
gresión.
1
Acrónimo de Akaike’s an Information Criterion
Capı́tulo 11: Modelos estadı́sticos en R 60
Debe tenerse en cuenta que la tabla ANOVA corresponde a una sucesión de modelos
ajustados. Las sumas de cuadrados que en ella aparecen corresponden a la disminución en
las sumas de cuadrados residuales como resultado de la inclusión de un término concreto
en un lugar concreto de la sucesión. Por tanto el orden de inclusión sólo será irrelevante en
experimentos ortogonales.
Para experimentos multiestrato el procedimiento consiste, en primer lugar, en proyectar
la respuesta sobre los estratos de error, una vez más en secuencia, y, después, en ajustar
la media del modelo a cada proyección. Para más detalles, consulte Chambers & Hastie
(1992).
Una alternativa más flexible a la tabla ANOVA completa es comparar dos o más modelos
directamente utilizando la función anova().
> anova([Link].1, [Link].2, ...)
El resultado es una tabla ANOVA que muestra las diferencias entre los modelos ajustados
cuando se ajustan en ese orden preciso. Los modelos ajustados objeto de comparación
constituyen por tanto una sucesión jerárquica. Este resultado no suministra información
distinta a la del caso predeterminado, pero facilita su comprensión y control.
2
Su acrónimo en inglés es ANOVA, de ANalysis Of Variance. En alguna literatura en español, el acrónimo
se sustituye por ANDEVA, de ANálisis DE VArianza.
Capı́tulo 11: Modelos estadı́sticos en R 61
conocida pero que puede variar con las observaciones, y µ es la media de y. Por
tanto, se supone que la distribución de y queda determinada por su media y tal vez un
parámetro de escala.
• La media, µ, es una función inversible del predictor lineal:
11.6.1 Familias
La clase de modelos lineales generalizados que pueden ser tratados en R incluye las
distribuciones de respuesta gaussian (normal), binomial, poisson, inverse gaussian (nor-
mal inversa) y gamma ası́ como los modelos de quasi-likelihood (cuasi-verosimilitud) cuya
distribución de respuesta no está explı́citamente definida. En este último caso debe especi-
ficarse la función de varianza como una función de la media, pero en el resto de casos esta
función está implı́cita en la distribución de respuesta.
Cada distribución de respuesta admite una variedad de funciones de enlace para conec-
tar la media con el predictor lineal. La tabla siguiente recoge las que están disponibles
automáticamente.
Aquı́ vamos a utilizar la segunda de estas convenciones, por lo que debemos añadir una
matriz a la hoja de datos:
> kalythos$Ymat <- cbind(kalythos$y, kalythos$n - kalythos$y)
Para ajustar los modelos utilizamos
> fmp <- glm(Ymat ∼ x, family = binomial(link=probit), data = kalythos)
> fml <- glm(Ymat ∼ x, family = binomial, data = kalythos)
Puesto que la función de enlace logit es la predeterminada, este parámetro puede omitirse
en la segunda expresión. Para ver los resultados de cada ajuste usaremos
> summary(fmp)
> summary(fml)
Ambos modelos se ajustan bien (demasiado bien). Para estimar LD50 podemos usar la
siguiente función:
> ld50 <- function(b) -b[1]/b[2]
> ldp <- ld50(coef(fmp)); ldl <- ld50(coef(fmp)); c(ldp, ldl)
y obtendremos los valores 43.663 y 43.601 respectivamente.
θ 1 z1
y= +e
z2 − θ 2
1
y= +e
β1 x1 + β2 x2
donde x1 = z2 /z1 , x2 = −1/x1 , β1 = 1/θ1 y β2 = θ2 /θ1 . Suponiendo que existe una hoja de
datos apropiada, bioquimica, podemos ajustar este modelo mediante
> nlfit <- glm(y ∼ x1 + x2 - 1,
family = quasi(link=inverse, variance=constant),
data = bioquimica)
Si desea mayor información, lea las ayudas correspondientes.
> x <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22, 0.56, 0.56,
1.10, 1.10)
> y <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200)
El modelo a ajustar es:
> fn <- function(p) sum((y - (p[1] * x)/(p[2] + x))^2)
Para realizar el ajuste necesitamos unos valores iniciales de los parámetros. Una forma de
encontrar unos valores iniciales apropiados es representar gráficamente los datos, conjeturar
unos valores de los parámetros, y dibujar sobre los datos la curva correspondiente a estos
valores.
> plot(x, y)
> xajustado <- seq(.02, 1.1, .05)
> yajustado <- 200*xajustado/(.1+xajustado)
> lines(spline(xajustado,yajustado))
Aunque se podrı́a tratar de encontrar unos valores mejores, los valores obtenidos, 200 y
0.1, parecen adecuados. Ya podemos realizar el ajuste:
> resultado <- nlm(fn,p=c(200,.1),hessian=TRUE)
Tras el ajuste, resultado$minimum contiene SSE, y resultado$estimates contiene los
estimadores por mı́nimos cuadrados de los parámetros. Para obtener los errores tı́picos
aproximados de los estimadores (SE) escribimos lo siguiente:
> sqrt(diag(2*resultado$minimum/(length(y) - 2) *
solve(resultado$hessian)))
El número 2 en dicha expresión representa el número de parámetros. Un intervalo de
confianza al 95% será: El estimador del parámetro ± 1.96 SE. Podemos representar el ajuste
en un nuevo gráfico:
> plot(x,y)
> xajustado <- seq(.02,1.1,.05)
> yajustado <- 212.68384222*xajustado/(0.06412146+xajustado)
> lines(spline(xajustado,yajustado))
La biblioteca nls contiene muchas más posibilidades para ajustar modelos no lineales por
mı́nimos cuadrados. El modelo que acabamos de ajustar es el modelo de Michaelis-Menten,
por tanto podemos usar
> df <- [Link](x=x, y=y)
> fit <- nls(y ∼ SSmicmen(x, Vm, K), df)
> fit
Nonlinear regression model
model: y ∼ SSmicmen(x, Vm, K)
data: df
Vm K
212.68370711 0.06412123
residual sum-of-squares: 1195.449
> summary(fit)
Parameters:
Estimate Std. Error t value Pr(>|t|)
Capı́tulo 11: Modelos estadı́sticos en R 67
12 Procedimientos gráficos
Las posibilidades gráficas son un componente de R muy importante y versátil. Es posi-
ble utilizarlas para mostrar una amplia variedad de gráficos estadı́sticos y también para
construir nuevos tipos de gráficos.
Los gráficos pueden usarse tanto en modo interactivo como no interactivo, pero en la
mayorı́a de los casos el modo interactivo es más productivo. Además al iniciar R en este
modo, se activa un dispositivo para mostrar gráficos. Aunque este paso es automático, es
útil conocer que la orden es X11(), aunque también puede usar windows() en Microsoft
Windows.
Las órdenes gráficas se dividen en tres grupos básicos:
• Alto nivel Son funciones que crean un nuevo gráfico, posiblemente con ejes, etiquetas,
tı́tulos, etc..
• Bajo nivel Son funciones que añaden información a un gráfico existente, tales como
puntos adicionales, lı́neas y etiquetas.
• Interactivas Son funciones que permiten interactuar con un gráfico, añadiendo o elimi-
nando información, utilizando un dispositivo apuntador, como un ratón.
Además, en R existe una lista de parámetros gráficos que pueden utilizarse para adaptar
los gráficos.
plot(hd)
plot(∼ expr)
plot(y ∼ expr)
Sean hd, una hoja de datos; y, un objeto cualquiera; y expr, una lista de
nombres de objetos separados por sı́mbolos ‘+’ (por ejemplo, a + b + c). Las
dos primeras formas producen diagramas de todas las parejas de variables de
la hoja de datos hd (en el primer caso) y de los objetos de la expresión expr
(en el segundo caso). La tercera forma realiza sendos gráficos de y sobre cada
objeto nombrado en la expresión expr (uno para cada objeto).
R posee dos funciones muy útiles para representar datos multivariantes. Si X es una
matriz numérica o una hoja de datos, la orden
> pairs(X)
produce una matriz de gráficos de dispersión para cada pareja de variables definidas por
las columnas de X, esto es, cada columna de X se representa frente a cada una de las demás
columnas, y los n(n−1) gráficos se presentan en una matriz de gráficos con escalas constantes
sobre las filas y columnas de la matriz.
Cuando se trabaja con tres o cuatro variables, la función coplot puede ser más apro-
piada. Si a y b son vectores numéricos y c es un vector numérico o un factor (todos de la
misma longitud) entonces la orden
> coplot(a ∼ b | c)
produce diagramas de dispersión de a sobre b para cada valor de c. Si c es un factor,
esto significa que a se representa sobre b para cada nivel de c. Si c es un vector numérico,
entonces se agrupa en intervalos y para cada intervalo se representa a sobre b para los valores
de c dentro del intervalo. El número y tamaño de los intervalos puede controlarse con el
argumento [Link] de la función coplot(). La función [Link]() también es
útil para seleccionar intervalos. Asimismo, es posible utilizar dos variables condicionantes
con una orden como
> coplot(a ∼ b | c + d)
que produce diagramas de a sobre b para cada intervalo de condicionamiento de c y d.
Las funciones coplot() y pairs() utilizan el argumento panel para personalizar el tipo
de gráfico que aparece en cada panel o recuadro. El valor predeterminado es points() para
producir un diagrama de dispersión, pero si se introducen otras funciones gráficas de nivel
bajo de los dos vectores x e y como valor de panel, se produce cualquier tipo de gráfico que
se desee. Una función útil en este contexto es [Link]().
Existen otras funciones gráficas de nivel alto que producen otros tipos de gráficos. Al-
gunas de ellas son las siguientes:
Capı́tulo 12: Procedimientos gráficos 71
qqnorm(x)
qqline(x)
qqplot(x, y)
Gráficos de comparación de distribuciones. El primero representa el vector x
sobre los valores esperados normales. El segundo le añade una recta que pasa
por los cuartiles de la distribución y de los datos. El tercero representa los
cuantiles de x sobre los de y para comparar sus distribuciones respectivas.
hist(x)
hist(x, nclass=n)
hist(x, breaks=b, ...)
Produce un histograma del vector numérico x. El número de clases se cal-
cula habitualmente de modo correcto, pero puede elegir uno con el argumento
nclass, o bien especificar los puntos de corte con el argumento breaks. Si está
presente el argumento probability=TRUE, se representan frecuencias relativas
en vez de absolutas.
dotplot(x, ...)
Construye un gráfico de puntos de x. En este tipo de gráficos, el eje y etiqueta
los datos de x y el eje x da su valor. Por ejemplo, permite una selección visual
sencilla de todos los elementos con valores dentro de un rango determinado.
image(x, y, z, ...)
contour(x, y, z, ...)
persp(x, y, z, ...)
Gráficos tridimensionales. image representa una retı́cula de rectángulos con
colores diferentes según el valor de z, contour representa curvas de nivel de z,
y persp representa una superficie tridimensional de z.
abline(a, b)
abline(h=y)
abline(v=x)
abline([Link])
Añade al gráfico actual, una recta de pendiente b y ordenada en el origen a.
La forma h=y representa una recta horizontal de altura y, y la forma v=x,
una similar, vertical. En el cuarto caso, [Link] puede ser una lista con un
componente coefficients de longitud 2 (como el resultado de una función de
ajuste de un modelo) que se interpretan como ordenada y pendiente, en ese
orden.
polygon(x, y, ...)
Añade al gráfico actual un polı́gono cuyos vértices son los elementos de (x,y);
(opcionalmente) sombreado con lı́neas, o relleno de color si el periférico lo ad-
mite.
legend(x, y, letrero, ...)
Añade al gráfico actual un letrero o leyenda, en la posición especificada. Los
caracteres para dibujar, los estilos de lı́neas, los colores, etc. están identificados
con los elementos del vector letrero. Debe darse al menos otro argumento
más, v, un vector de la misma longitud que letrero, con los correspondientes
valores de dibujo, como sigue:
legend( , fill=v)
Colores para rellenar
legend( , col=v)
Colores de puntos y lı́neas
legend( , lty=v)
Tipos de lı́nea
legend( , lwd=v)
Anchura de lı́nea
legend( , pch=v)
Caracteres para dibujar (vector de caracteres)
title(main, sub)
Añade un tı́tulo, main, en la parte superior del gráfico actual, de tamaño grande,
y un subtı́tulo, sub, en la parte inferior, de tamaño menor.
axis(side, ...)
Añade al gráfico actual un eje en el lado indicado por el primer argumento (de
1 a 4, en el sentido de las agujas del reloj, siendo el 1 la parte inferior). Otros
argumentos controlan la posición de los ejes, dentro o fuera del gráfico, las
marcas y las etiquetas. Es útil para añadir ejes tras utilizar la función plot()
con el argumento axes=FALSE.
Las funciones gráficas de nivel bajo necesitan normalmente alguna información de
posición, como las coordenadas x e y, para determinar dónde colocar los nuevos elementos.
Las coordenadas se dan en términos de coordenadas de usuario, las cuales están definidas
Capı́tulo 12: Procedimientos gráficos 74
por las funciones gráficas de alto nivel previas y se toman en función de los datos
suministrados a estas funciones.
Los dos argumentos, x e y, pueden sustituirse por un solo argumento de clase lista con dos
componentes llamados x e y, o por una matriz con dos columnas. De este modo, funciones
como locator(), que vemos a continuación, pueden usarse para especificar interactivamente
posiciones en un gráfico.
nombre=valor
Descripción del efecto del parámetro. nombre es el nombre del parámetro, esto
es, el nombre del argumento que debe usar en la función par() o cualquier
función gráfica. valor es un valor tı́pico del parámetro.
3.0
Plot region
1.5
0.0
y
mai[2]
−1.5
−3.0
mai[1] x
Margin
mar=c(4, 2, 2, 1)
Similar a mai, pero medido en lı́neas de texto.
Los parámetros mar y mai están relacionados en el sentido de que un cambio en uno se
refleja en el otro. Los valores predeterminados son a menudo demasiado grandes, el margen
derecho se necesita raramente, igual que el superior si no se incluye tı́tulo. Los márgenes
inferior e izquierdo sólo necesitan el tamaño preciso para incluir las etiquetas de ejes y de
división. Además el valor predeterminado no tiene en cuenta la superficie del dispositivo
gráfico. Ası́, si utiliza el dispositivo postscript() con el argumento height=4 obtendrá
un gráfico en que la mitad del mismo son márgenes, salvo que explı́citamente cambie mar
o mai. Cuando hay figuras múltiples, como veremos después, los márgenes se reducen a la
mitad, aunque suele ser insuficiente cuando varias figuras comparten la misma página.
R permite la creación de una matriz de n × m figuras en una sola página. Cada figura
tiene sus propios márgenes, y la matriz de figuras puede estar opcionalmente rodeada de
un margen exterior , tal como se muestra en la siguiente figura:
Capı́tulo 12: Procedimientos gráficos 80
−−−−−−−−−−−−−−−
−−−−−−−−−−−−−−−
−−−−−−−−−−−−−−− oma[3]
−−−−−−−−−−−−−−−
−−−−−−−−−−−−−−−
omi[4]
mfg=c(3,2,3,2)
omi[1]
mfrow=c(3,2)
Los parámetros gráficos relacionados con las figuras múltiples son los siguientes:
mfcol=c(3, 2)
mfrow=c(2, 4)
Definen el tamaño de la matriz de figuras múltiples. En ambos, el primer valor
es el número de filas, el segundo el de columnas. La diferencia entre los dos, es
que con el primero, mfcol, la matriz se rellena por columnas, en tanto que con
el segundo, mfrow, lo hace por filas. La distribución en la figura del ejemplo se
ha creado con mfrow=c(3,2) y se ha representado la página en el momento en
que se han realizado los cuatro primeros gráficos.
mfg=c(2, 2, 3, 2)
Definen la posición de la figura actual dentro de la matriz de figuras múltiples.
Los dos primeros valores indican la fila y columna de la figura actual, en tanto
que los dos últimos son el número de filas y columnas de la matriz de figuras
múltiples. Estos parámetros se utilizan para seleccionar cada una de las dife-
rentes figuras de la matriz. Incluso, los dos últimos valores pueden ser distintos
de los verdaderos valores, para poder obtener figuras de tamaños distintos en
la misma página.
fig=c(4, 9, 1, 4)/10
Definen la posición de la figura actual en la página. Los valores son las posiciones
de los bordes izquierdo, derecho, inferior y superior respectivamente, medidos
como la proporción de página desde la esquina inferior izquierda. El ejemplo
corresponderı́a a una figura en la parte inferior derecha de la página. Este
parámetro permite colocar una figura en cualquier lugar de la página.
oma=c(2, 0, 3, 0)
omi=c(0, 0, 0.8, 0)
Definen el tamaño de los márgenes exteriores. De modo similar a mar y mai, el
primero está expresado en lı́neas de texto y el segundo en pulgadas, correspon-
den a los márgenes inferior, izquierdo, superior y derecho respectivamente.
Los márgenes exteriores son particularmente útiles para ajustar convenientemente los
tı́tulos, etc. Puede añadir texto en estos márgenes con la función mtext() sin más que
Capı́tulo 12: Procedimientos gráficos 81
consecuencia de la compatibilidad con S e indica que la salida está constituida por una sola
página. Por tanto para crear un gráfico que se pueda incluir sin problema en cualquier
procesador, deberá utilizar una orden análoga a la siguiente:
> postscript("[Link]", horizontal=FALSE, onefile=FALSE,
height=8, width=6, pointsize=10)
windows()
Abre una ventana en Microsoft Windows
postscript()
pictex()
... Cada llamada a una función de controlador de dispositivo abre un nuevo dis-
positivo gráfico y, por tanto, añade un elemento a la lista de dispositivos, al
tiempo que este dispositivo pasa a ser el dispositivo actual al que se enviarán
los resultados gráficos. (En algunas plataformas es posible que existan otros
dispositivos disponibles.)
[Link]()
Devuelve el número y nombre de todos los dispositivos activos. El dispositivo
de la posición 1 es siempre el dispositivo nulo que no acepta ninguna orden
gráfica.
[Link]()
[Link]()
Devuelve el número y nombre del dispositivo gráfico siguiente o anterior, re-
spectivamente, al dispositivo actual.
[Link](which=k)
Puede usarse para hacer que el dispositivo que ocupa la posición k en la lista
de dispositivos sea el actual. Devuelve el número y nombre del dispositivo.
[Link](k)
Cierra el dispositivo gráfico que ocupa la posición k de la lista de dispositivos.
Para algunos dispositivos, como los postscript, finalizará el gráfico, bien im-
primiéndolo inmediatamente, bien completando la escritura en el archivo, de-
pendiendo de cómo se ha iniciado el dispositivo.
Capı́tulo 12: Procedimientos gráficos 83
attach(mm)
Conecta la hoja de datos en la posición predeterminada: la 2.
plot(Expt, Speed, main="Velocidad de la luz", xlab="No. Experimento")
Compara los cinco experimentos mediante diagramas de cajas.
fm <- aov(Speed ∼ Run + Expt, data=mm)
summary(fm)
Analiza los datos como un diseño en bloques aleatorizados, tomando las ‘series’
y los ‘experimentos’ como factores.
fm0 <- update(fm, . ∼ . - Run)
anova(fm0,fm)
Ajusta el submodelo eliminando las ‘series’, y lo compara utilizando un análisis
de la varianza.
detach()
rm(fm, fm0)
Desconecta la hoja de datos y elimina los objetos creados, antes de seguir ade-
lante.
Consideraremos ahora nuevas posibilidades gráficas.
x <- seq(-pi, pi, len=50)
y <- x x es un vector con 50 valores equiespaciados en el intervalo −π ≤ x ≤ π. El
vector y es idéntico a x.
f <- outer(x, y, function(x, y) cos(y)/(1 + x^2))
f es una matriz cuadrada, cuyas filas y columnas están indexadas por x e y
respectivamente, formada por los valores de la función cos(y)/(1 + x2 ).
oldpar <- par([Link] = TRUE)
par(pty="s")
Almacena los parámetros gráficos y modifica el parámetro pty (zona de dibujo)
para que valga “s” (cuadrado).
contour(x, y, f)
contour(x, y, f, nlevels=15, add=TRUE)
Dibuja un mapa de curvas de nivel de f ; y después le añade más lı́neas para
obtener más detalle.
fa <- (f-t(f))/2
fa es la “parte asimétrica” de f . (t(f) es la traspuesta de f).
contour(x, y, fa, nint=15)
Dibuja un mapa de curvas de nivel,. . .
par(oldpar)
. . . y recupera los parámetros gráficos originales.
image(x, y, f)
image(x, y, fa)
Dibuja dos gráficos de densidad
Apendice A: Primera sesión con R 87
Apendice B Ejecución de R
--help
-h Muestra un pequeño mensaje de ayuda y termina correctamente.
--version
Muestra la información de la versión y termina correctamente.
RHOME Muestra la trayectoria al “directorio inicial” de R y termina correctamente.
Salvo la página de ayuda en UNIX y el archivo de ejecución, la instalación de
R pone todos los archivos (ejecutables, bibliotecas, etc.) en este directorio.
--save
--no-save
Indica si debe salvar la imagen del entorno al terminar la sesión. En modo in-
teractivo, si no ha especificado ninguna, le preguntará. En modo no interactivo
es obligatorio usar una.
--no-environ
No lee ningún archivo para dar valor a las variables de entorno.
--no-site-file
No lee el perfil de inicio global al iniciar el programa.
--no-init-file
No lee el perfil de inicio de usuario al iniciar el programa.
--restore
--no-restore
Indica si la imagen salvada (archivo ‘.Rdata’ en el directorio en que R se haya
iniciado) debe ser recuperada o no. El valor predeterminado es recuperarla.
--vanilla
Combina las opciones ‘--no-save’, ‘--no-environ’ ‘--no-site-file’,
‘--no-init-file’, y ‘--no-restore’.
--no-readline
Desactiva la edición de órdenes a través de readline. Esta opción suele utilizarse
cuando se ejecuta R desde Emacs utilizando la biblioteca ESS (“Emacs Speaks
Statistics”). Véase Apendice C [El editor de ordenes], página 93, para ampliar
esta información.
--vsize=N
Indica la cantidad de memoria utilizada para objetos de tamaño variable. N
debe ser un entero, en cuyo caso la unidad de medida es el octeto o byte, o
un entero seguido de una letra que indica la unidad de medida en octetos: ‘M’,
‘K’, o ‘k’, que corresponden respectivamente a ‘Mega’ (2^20), ‘Kilo’ (2^10) (de
ordenador), o ‘kilo’decimal (1000).
--nsize=N
Indica la cantidad de memoria utilizada para objetos de tamaño fijo. Las con-
sideraciones realizadas en el apartado anterior son válidas para N.
--quiet
--silent
-q No muestra ni el mensaje de copyright ni los mensajes de inicio.
Apendice B: Ejecución de R 90
--slave Ejecuta R con el mı́nimo de salidas posible. Esta opción se utiliza con programas
que se sirven de R para realizar cálculos para ellos mismos.
--verbose
Muestra el máximo de salidas posible. Además modifica la opción verbose a
TRUE. R utiliza esta opción para controlar si debe mostrar mensajes de diag-
nóstico.
--debugger=depurador
-d depurador
Ejecuta R desde el programa de depuración (debugger) depurador. En este
caso, si existen otras opciones, se descartan. Cualquier otra opción debe darse
al iniciar R desde el programa de depuración.
--gui=tipo
Utiliza tipo como interfaz gráfico (advierta que esto también incluye los gráficos
interactivos). Los valores posibles de tipo son X11 (predeterminado) y GNOME,
supuesto que GNOME esté disponible.
La entrada y la salida pueden redirigirse de la manera habitual, utilizando ‘<’ and ‘>’.
R CMD permite utilizar varias herramientas que son útiles en conjunción con R, pero que
no están diseñadas para usarlas “directamente”. La forma general es
R CMD orden argumentos
donde orden es el nombre de la herramienta y argumentos son los argumentos que se pasan
a la misma.
Las herramientas disponibles son:
BATCH Ejecuta R en modo no interactivo.
COMPILE Compila archivos para uso con R.
SHLIB Construye bibliotecas compartidas del sistema operativo para carga dinámica.
INSTALL Instala bibliotecas añadidas.
REMOVE Elimina bibliotecas añadidas.
build Construye bibliotecas añadidas.
check Comprueba bibliotecas añadidas.
Rdconv Convierte desde formato Rd a otros formatos, incluyendo html, Nroff, LaTEX,
texto ASCII sin formato, y formato de documentación S.
Rd2dvi Convierte desde formato Rd a DVI/PDF.
Rd2txt Convierte desde formato Rd a texto con formato.
Rdindex Extrae información de ı́ndices de archivos Rd.
Sd2Rd Convierte desde formato de documentación S a formato Rd.
Las cinco primeras herramientas (BATCH, COMPILE, SHLIB, INSTALL, y REMOVE) pueden
ejecutarse “directamente” sin la opción CMD, esto es, en la forma R orden argumentos.
Utilice la orden
R CMD herramienta --help
para obtener información sobre cada una de las herramientas descritas.
Apendice B: Ejecución de R 91
Las opciones que pueden darse al ejecutar R en Microsoft Windows son las siguientes:
--version
Muestra la información de la versión y termina correctamente.
--mdi
--sdi
--no-mdi Inidica si Rgui se comportará como un programa MDI (predeterminado), donde
cada nueva ventana está contenida dentro de la ventana principal, o como un
programa SDI, donde cada ventana aparece de modo independiente en el es-
critorio.
--save
--no-save
Indica si debe salvar la imagen del entorno al terminar la sesión. En modo in-
teractivo, si no ha especificado ninguna, le preguntará. En modo no interactivo
es obligatorio usar una.
--restore
--no-restore
Indica si la imagen salvada (archivo ‘.Rdata’ en el directorio en que R se haya
iniciado) debe ser recuperada o no. El valor predeterminado es recuperarla.
--no-site-file
No lee el perfil de inicio global al iniciar el programa.
--no-init-file
No lee el archivo ‘.Rprofile’ del directorio del usuario (perfil de inicio de
usuario) al iniciar el programa.
--no-environ
No lee el archivo ‘.Renviron’.
--vanilla
Combina las opciones --no-save, --no-restore, --no-site-file,
--no-init-file y --no-environ.
-q
--quiet
--silent No muestra el mensaje de inicio.
--slave Ejecuta R con el mı́nimo de salidas posible.
--verbose
Muestra el máximo de salidas posible.
--ess Prepara Rterm para uso en modo R-inferior en ESS.
Apendice C: El editor de órdenes 93
C.1 Preliminares
Si la biblioteca de GNU, ‘readline’, está disponible cuando se compila R en UNIX, se
puede utilizar un editor interno de lı́neas de órdenes que permite recuperar, editar y volver
a ejecutar las órdenes utilizadas previamente.
Este editor puede desactivarse (lo que permite utilizar ESS1 ) dando la opción de inicio
‘--no-readline’.
La versión de Microsoft Windows dispone de un editor más sencillo. Vea la opción
‘Console’ en el menú ‘Help’.
Cuando use R con las posibilidades de readline, estarán disponibles las opciones que
posteriormente se indican.
Tenga en cuenta las siguientes convenciones tipográficas usuales:
Muchas de las órdenes utilizan caracteres caracteres Control y Meta. Los caracteres
Control, tales como Control-m, se obtienen pulsando la tecla hCTRLi y, sin soltarla, pulsando
la tecla hmi, y en la tabla los escribiremos como C-m. Los caracteres Meta, tales como Meta-
b, se obtienen pulsando la tecla hMETAi y, después de soltarla, pulsando la tecla hbi, y en la
tabla los escribiremos como M-b. Si su teclado no tiene la tecla hMETAi, puede obtenerlos
mediante una secuencia de dos caracteres que comienza con ESC. Esto es, para obtener M-b,
deberá escribir hESCihbi. Estas secuencias ESC también puede utilizarlas aunque su teclado sı́
disponga de la tecla hMETAi. Debe tener en cuenta que en los caracteres Meta se distingue
entre mayúsculas y minúsculas, por lo que puede que sea distinto el resultado si se pulsa
M-b o M-B.
1
Corresponde al acrónimo del editor de textos, ‘Emacs Speaks Statistics’; vea la dirección
[Link]
Apendice C: El editor de órdenes 94
Edición
text Inserta texto en el cursor.
C-f text Añade texto tras el cursor.
hDELi Borra el carácter a la izquierda del cursor.
C-d Borra el carácter bajo el cursor.
M-d Borra el resto de la palabra bajo el cursor, y la guarda.
C-k Borra el resto de la lı́nea desde el cursor, y lo guarda.
C-y Inserta el último texto guardado.
C-t Intercambia el carácter bajo el cursor con el siguiente.
M-l Cambia el resto de la palabra a minúsculas.
M-c Cambia el resto de la palabra a mayúsculas.
hRETi Vuelve a ejecutar la lı́nea.
Al pulsar hRETi, se termina la edición de la lı́nea.
Apendice D: Índice de funciones y variables 95
! +
! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
+ ............................................ 8
!= . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
% >
%*% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
%o% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 > . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
>= . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
&
& . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
&& . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 ^
^ ............................................ 8
*
* ............................................ 8
<
-
< . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
- ............................................ 8
<= . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
<<- . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
.First . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
.Last . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
A
abline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
/ ace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
/ ............................................ 8 add1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
anova . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59, 60
aov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
:
aperm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
: ............................................ 9
array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
= [Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
== . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 attach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
attr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
? attributes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
? ............................................ 4 avas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
|
| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
|| . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 B
boxplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
~ break . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
~ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 bruto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Apendice D: Índice de funciones y variables 96
C G
c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7, 10, 27, 29 glm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
cbind . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
coef . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
H
coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 help . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
contour . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 hist . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38, 71
Contrastes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
coplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
cos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
I
crossprod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22, 24 identify. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
cut . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 if . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
ifelse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
D [Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 K
density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 [Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
detach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
[Link]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
[Link]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 L
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 legend . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
[Link]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8, 14
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
deviance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
diag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 list . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
dim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 lm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
dotplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 lme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
drop1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 locator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
loess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67, 68
log . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
lqs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
E lsfit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
ecdf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
eigen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
else . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 M
Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 mars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
exp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 max . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
mean . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
min . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
F mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
F ............................................ 9
factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 N
FALSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 NA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
fivenum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 NaN . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
for . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 ncol . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 next . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 nlm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65, 66, 67
Apendice D: Índice de funciones y variables 97
nlme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 S
nrow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 scan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
seq . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
O [Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
order . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 sin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
ordered . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 sink . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
sort . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
outer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
source . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
split . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
sqrt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
P stem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
pairs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59, 61
par . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 sum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
paste . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38, 59
persp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 svd . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
pictex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59, 69 T
pmax . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
pmin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
T ............................................ 9
points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
polygon . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 table . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22, 27
postscript . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 tan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
predict . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 tapply . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
print . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 text . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
prod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 title . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
TRUE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
Q
qqline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39, 71 U
qqnorm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39, 71 unclass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
qqplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
qr . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
V
R var . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
range . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 [Link]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
vector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
rbind . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
[Link]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
[Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 W
rep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 while . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
repeat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 [Link] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
resid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
residuals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
rlm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 X
rm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 X11 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
Apendice E: Índice de conceptos 98
A G
Acceso a datos internos . . . . . . . . . . . . . . . . . . . . . . 35 Gráficos cuantil-cuantil . . . . . . . . . . . . . . . . . . . . . . . 39
Actualización de modelos ajustados . . . . . . . . . . . 61 Gráficos dinámicos . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
Ajuste por mı́nimos cuadrados . . . . . . . . . . . . . . . . 25
Ámbito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
Análisis de varianza . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Argumentos con nombre. . . . . . . . . . . . . . . . . . . . . . 47
H
Asignación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Hojas de datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
Atributos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Autovalores y autovectores . . . . . . . . . . . . . . . . . . . 25
I
B Indexación de vectores . . . . . . . . . . . . . . . . . . . . . . . 11
Bibliotecas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Indexación de y mediante variables indexadas . . 20
C K
Ciclos y ejecución condicional . . . . . . . . . . . . . . . . . 44
Kolmogorov-Smirnov, contraste de . . . . . . . . . . . . 40
Clases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15, 54
Cómo definir un operador binario . . . . . . . . . . . . . 47
Cómo importar datos . . . . . . . . . . . . . . . . . . . . . . . . 36
Concatenación de listas . . . . . . . . . . . . . . . . . . . . . . 29 L
Contrastes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 Lectura de datos desde un archivo . . . . . . . . . . . . 33
Contrastes de una y de dos muestras . . . . . . . . . . 41 Listas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
Controladores de dispositivos gráficos . . . . . . . . . 81
D M
Descomposición en valores singulares . . . . . . . . . . 25 Máxima verosimilitud . . . . . . . . . . . . . . . . . . . . . . . . 67
Descomposición QR . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Mı́nimos cuadrados no lineales . . . . . . . . . . . . . . . . 65
Determinantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 Modelos aditivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Diagrama de cajas . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 Modelos basados en árboles . . . . . . . . . . . . . . . . . . . 68
Distribuciones de probabilidad . . . . . . . . . . . . . . . . 37 Modelos estadı́sticos . . . . . . . . . . . . . . . . . . . . . . . . . 55
Modelos lineales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
Modelos lineales generalizados . . . . . . . . . . . . . . . . 61
E Modelos mezclados. . . . . . . . . . . . . . . . . . . . . . . . . . . 67
Eliminar objetos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Escritura de funciones . . . . . . . . . . . . . . . . . . . . . . . . 46
Espacio de trabajo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
Estimación de la densidad . . . . . . . . . . . . . . . . . . . . 38
O
Expresiones agrupadas . . . . . . . . . . . . . . . . . . . . . . . 44 Objetos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
Órdenes de control . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
Orientación a objetos . . . . . . . . . . . . . . . . . . . . . . . . 54
F
Factores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17, 58
Factores ordinales . . . . . . . . . . . . . . . . . . . . . . . . 17, 58 P
Familias. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
Fórmulas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Parámetros gráficos . . . . . . . . . . . . . . . . . . . . . . . . . . 76
Función de distribución empı́rica . . . . . . . . . . . . . . 39 Personalización del entorno . . . . . . . . . . . . . . . . . . . 53
Funciones genéricas . . . . . . . . . . . . . . . . . . . . . . . . . . 54 Producto exterior de variables indexadas . . . . . . 23
Funciones y operadores aritméticos . . . . . . . . . . . . . 8 Producto matricial . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
Apendice E: Índice de conceptos 99
R Tabulación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
Redirección de la entrada y la salida . . . . . . . . . . . 5 Traspuesta generalizada de una variable indexada
Regla de reciclado . . . . . . . . . . . . . . . . . . . . . . . . . 8, 22 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
Regresión con aproximación local . . . . . . . . . . . . . 67 Trayectoria de búsqueda. . . . . . . . . . . . . . . . . . . . . . 31
Regresión robusta. . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
V
S Valores faltantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
Shapiro-Wilk, contraste de . . . . . . . . . . . . . . . . . . . 40 Valores predeterminados . . . . . . . . . . . . . . . . . . . . . 47
Student, contraste t de . . . . . . . . . . . . . . . . . . . . . . . 41 Vectores de caracteres . . . . . . . . . . . . . . . . . . . . . . . . 10
Sucesiones regulares . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
W
T Wilcoxon, contraste de . . . . . . . . . . . . . . . . . . . . . . . 42
Apendice F: Referencias 100
Apendice F Referencias
D. M. Bates and D. G. Watts (1988), Nonlinear Regression Analysis and Its Applications.
John Wiley & Sons, New York.
Richard A. Becker, John M. Chambers and Allan R. Wilks (1988), The New S Language.
Chapman & Hall, New York. This book is often called the “Blue Book ”.
John M. Chambers and Trevor J. Hastie eds. (1992), Statistical Models in S. Chapman
& Hall, New York. This is also called the “White Book ”.
Annette J. Dobson (1990), An Introduction to Generalized Linear Models, Chapman and
Hall, London.
Peter McCullagh and John A. Nelder (1989), Generalized Linear Models. Second edition,
Chapman and Hall, London.
John A. Rice (1995), Mathematical Statistics and Data Analysis. Second edition.
Duxbury Press, Belmont, CA.
S. D. Silvey (1970), Statistical Inference. Penguin, London.
Econometría básica con Python .
Esta obra está bajo la licencia Creative Commons Atribución-NoComercial 4.0 Internacional
Prefacio
Esta obra constituye un manual para la realización de análisis econométrico básico con el uso del lenguaje de
programación Python; no comprende, como parte de su contenido, un cuerpo teórico sólido que pueda tomarse
como base adecuada para la enseñanza de la Econometría y tan solo se limita a ser un referente práctico de la
aplicación de las técnicas y modelos propios de la materia. Dado el enfoque aplicado de este trabajo, el lector
no encontrará en éste un texto que le sirva de guía teórica y le permita fortalecer sus conocimientos sobre la
Econometría, pues no se trata del objetivo aquí planteado; para tal propósito, se sugiere recurrir a los manuales
especializados generalmente empleados como material de referencia en los cursos de pregrado en Economía.
El contenido que se presenta en esta obra está especialmente dirigido a estudiantes de Economía, Estadís-
tica y carreras afines, a los profesionales cuyas labores se vean asociadas al ámbito de la Econometría y a toda
persona, con cierto conocimiento previo, que sienta afinidad por la materia y esté interesada en mejorar sus ha-
bilidades en el uso de herramientas informáticas útiles en el desarrollo de ejercicios de econometría aplicada.
Dado el mínimo nivel de profundidad con el que se abordan los fundamentos teóricos sobre los que se estruc-
tura el análisis econométrico en esta obra, y su marcado enfoque práctico, se espera que el lector posea unos
buenos conocimientos previos (a nivel introductorio) sobre la materia, de modo que no le resulte complejo
seguir los ejercicios expuestos y pueda alcanzar sus objetivos de aprendizaje.
Además de un conocimiento básico previo sobre Econometría, se espera que el lector esté familiarizado
con el uso de herramientas informáticas, de forma que comprenda, sin mayores problemas, las indicaciones
dadas para el desarrollo de los ejercicios planteados y las instrucciones brindadas para la instalación y uso
de los programas empleados. Asimismo, resulta necesario que el lector posea cierto conocimiento del idioma
inglés, en cuanto la mayoría de funciones, métodos y atributos empleados son representados por medio de
palabras en inglés, las cuales, dado el muy alto nivel del lenguaje de programación que se empleará, resultan
muy informativas y facilitan enormemente la comprensión de las tareas que se llevan a cabo.
Deseo resaltar mi interés en que este trabajo se constituya como referente no teórico y facilite el aprendiza-
je de la Econometría dado el tratamiento práctico que pretende lograr; asimismo, que permita el acercamiento
de toda persona interesada en la materia, con cierto nivel de conocimiento, sin importar el ámbito de su for-
mación, y que promueva el aprovechamiento de herramientas tecnológicas por parte de todos. Aunque cierto
conocimiento o experiencia previa con lenguajes de programación son bienvenidos, y resultarán muy útiles para
la comprensión de la obra por parte del lector, estos no son indispensables pues, con el tratamiento dado a la
información en esta obra, pretendo que las exposiciones presentadas sean lo más ilustrativas y claras posibles,
de modo que todo lector se sienta cómodo con el estudio de este trabajo y alcance sus metas de aprendizaje
satisfactoriamente.
El Autor,
Fabián Alejandro Triana Alarcón
Índice general
4. Verificación de supuestos 39
Supuestos sobre la estructura del modelo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Preparación del entorno . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Importación de los datos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
Número de observaciones mayor a número de parámetros . . . . . . . . . . . . . . . . . . . . 41
Variación en las variables explicativas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
No multicolinealidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
No sesgo de especificación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Supuestos sobre el término de error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Valor medio igual a cero . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
Homoscedasticidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
No autocorrelación . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
Normalidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
Referencias 111
Capítulo 1
A lo largo de este trabajo se pretende desarrollar ejercicios prácticos, con propósitos ilustrativos, que sirvan
de referente para el uso del lenguaje de programación Python como herramienta útil para el análisis economé-
trico y que permitan al lector familiarizarse con el mismo y fortalecer sus habilidades en el manejo de éste.
Esta obra se basa en el uso de Python 3, con distribución Anaconda y trabajando sobre Jupyter Notebook;
aunque estos nombres resulten extraños para el lector, especialmente el novato, es adecuado señalar que éste
no debería preocuparse, ya que se trata de conceptos sin mayor complejidad (para nuestros propósitos) y que,
a continuación, se pretenden aclarar, de modo que sienta entera confianza en el estudio de esta obra.
Python
Python es un lenguaje de programación orientado a objetos, interpretado, de alto nivel y con tipado di-
námico (Python Software Foundation,s.f), creado por Guido van Rossum, un científico de la computación y
programador holandés, y que apareció en su versión inicial en 1991. Python es de muy alto nivel, lo que sugiere
que su sintaxis resulta relativamente fácil de comprender y lo que se refleja en el hecho de que su proceso de
aprendizaje es más bien sencillo en comparación con el que corresponde a otros lenguajes de programación.
Se ha señalado que trabajaremos con Python 3, sin embargo, no se ha informado al lector sobre la diferencia
entre Python, a secas, y Python 3. Esto no debería causarle intranquilidad, pues la distinción es bien simple:
Python es el lenguaje de programación, mientras que Python 3 es una versión de dicho lenguaje (el cual también
está disponible en la versión Python 2). En este trabajo, se hará uso de Python 3 ya que es la versión que
evidencia una tendencia creciente a la adopción generalizada y se caracteriza por un desarrollo activo.
El uso que se hace de Python a lo largo de esta obra es aplicado y específico, es decir, se emplea con propósitos
muy claramente delineados, ejecutando tareas útiles para abordar el análisis econométrico que se pretende
llevar a cabo; en ningún momento este trabajo trata de ser un manual de programación o algo que se le parezca,
por lo cual el lector no debe considerar que con el estudio de esta obra aprenderá a programar, pues tal objetivo
excede, y por mucho, el alcance de esta guía. Al lector interesado en la programación como tal, se sugiere
consultar otro tipo de recursos que le resulten útiles a tal propósito, en cuanto en esta obra se recurre a Python
únicamente como herramienta y no como el tema central a abordar.
Tal y como apenas se ha señalado, Python es una herramienta que utilizaremos para desarrollar análisis
econométrico; sin embargo, puede que tal planteamiento no resulte del todo claro al lector, quien es natural
que se cuestione ¿cómo se hará uso de Python para llevar a cabo análisis econométrico? La respuesta a tal
interrogante no implica mayor complejidad: tan solo se hará uso de las funciones, métodos y atributos que
nos ofrecen Python y sus librerías especializadas (las cuales abordaremos posteriormente) para manejar los
1
datos, extraer la información requerida, construir los modelos, realizar las estimaciones y ejecutar las pruebas
estadísticas necesarias. Aunque tales tareas sí comprenden operaciones con cierto nivel de dificultad e implican
cálculos elaborados, el usuario no tendrá que realizar un esfuerzo extraordinario, de hecho, ni siquiera tendrá
que definir una función de dificultad significativa por su propia cuenta o escribir código “complejo”, en cuanto
Python y sus librerías especializadas brindarán los instrumentos que requeriremos a lo largo de nuestra labor.
Anaconda
Aunque es posible instalar Python usando la distribución oficial del sitio web de Python Software Founda-
tion, [Link] no se recurrirá a tal medio, en cuanto se hará uso de la distribución Anaconda.
La explicación para tal decisión se encuentra en el hecho de que la distribución Anaconda instala automáti-
camente múltiples librerías (muy útiles para las labores que se adelantarán), incluye un gestor de paquetes
propio, conda, el sistema de gestión de paquetes estándar de Python (pip) y, además, nos da acceso a Anacon-
da Navigator, una interfaz gráfica con la que podemos ingresar a Jupyter Notebook, que es donde, en último
término, escribiremos el código Python con el que llevaremos a cabo nuestro análisis econométrico.
Anaconda Distribution es un producto de la empresa estadounidense Anaconda, Inc., (anteriormente cono-
cida como Continuum Analytics) y consiste en una distribución Python gratuita y de código abierto ampliamente
utilizada en el campo de la ciencia de datos y machine learning. Para mayor información, se sugiere consultar el
sitio web de la documentación oficial [Link] aunque, para nuestros propósi-
tos, esto no resultará ni siquiera necesario, en cuanto aquí se darán las instrucciones requeridas por el lector
para seguir correctamente el desarrollo de los temas tratados.
Jupyter Notebook
Jupyter Notebook es una aplicación web de código abierto que permite crear y compartir documentos que
contienen código, gráficos, ecuaciones y texto (Project Jupyter, 2019). Soporta varios lenguajes de programa-
ción, entre ellos Julia, Python y R, cuenta con una interfaz atractiva, de fácil manejo y muy amigable con el
usuario, y hace parte de Project Jupyter, una iniciativa surgida en 2014 a partir de IPython y dirigida por Fernan-
do Pérez, físico colombiano de la Universidad de Antioquia, profesor asistente del Departamento de Estadística
de UC Berkeley e investigador del Berkeley Institute for Data Science.
La elección de Jupyter Notebook para ser la aplicación en la que se desarrollará nuestro análisis economé-
trico radica en lo agradable de su interfaz y su manejo extremadamente simple e intuitivo, lo que nos permitirá
escribir código, visualizar resultados, elaborar gráficos y realizar comentarios en un mismo documento, enri-
queciendo vastamente nuestro análisis y facilitando la comprensión del contenido creado (de hecho, este libro
ha sido escrito en Jupyter Notebook, sobre LATEXy markdown, con edición posterior en Overleaf). Para mayor
información sobre Project Jupyter, se invita a consultar el sitio web oficial del proyecto: [Link]
2
Lo primero que debemos hacer es contar con la distribución Anaconda. En caso de no tenerla instalada
previamente (lo que es muy probable, especialmente para un usuario no familiarizado con el tema), debemos
acceder al sitio web oficial [Link] Al ingresar a tal dirección, nos encontra-
remos con una página como la siguiente:
Podemos observar que en la parte inferior se encuentra una sección en la que está escrito Windows | macOS
| Linux. En este punto, debemos hacer clic sobre el nombre correspondiente al sistema operativo de nuestro
computador; esto, con el propósito de acceder al Instalador Anaconda para nuestro sistema operativo. En es-
te caso particular, seleccionaremos Windows, sin embargo, el lector tiene que ser cuidadoso, en cuanto debe
seleccionar el sistema operativo de su computador (el cual, puede no ser Windows). Al seleccionar el sistema
operativo correspondiente, podremos visualizar lo siguiente:
Observamos que hay dos instaladores: uno para Python 3 (versión 3.71 específicamente) y otro para Python
2 (versión 2.72 específicamente); el que debemos seleccionar, como se ha señalado con anterioridad, es el de
1 La versión más reciente al momento de publicación de esta obra.
2 La versión más reciente al momento de publicación de esta obra.
3
Python 3. En cuanto al Graphical Installer específico, éste depende del procesador de nuestro computador: si
se trata de un procesador de 64 bits, naturalmente descargaremos el 64-Bit Graphical Installer, y si se trata de
un procesador de 32 bits, nuestra elección será el 32-Bit Graphical Installer. Para el lector que desconozca el
sistema de su computador, éste puede ser consultado, en Windows, en las propiedades del equipo o en:
Allí, encontrará, entre otros datos, la información correspondiente al sistema; información similar a la que
se presenta a continuación:
Una vez seleccionado el Graphical Installer correspondiente, solo se debe guardar el archivo y, una vez descarga-
do por completo, se debe ejecutar. Este proceso es muy simple (pues tan solo consiste en seguir instrucciones
muy concisas), por lo cual no se expondrá paso a paso, en cuanto hacerlo resultaría superfluo. Una vez des-
cargado e instalado el programa, accederemos a Anaconda Navigator; para tal propósito, recurriremos a una
búsqueda desde el botón de inicio de Windows, tal y como se señala a continuación:
4
Haremos clic sobre Anaconda Navigator (Aplicación de escritorio) y esperaremos que sea cargada por com-
pleto (esto puede llegar a tomar un poco de tiempo, por lo que el usuario no debe inquietarse). Recibiremos
entonces el siguiente mensaje:
5
Sobre este mensaje, simplemente debemos pulsar en Ok o en Ok, and don’t show again. Una vez rea-
lizado este paso, tendremos acceso a Anaconda Navigator, donde encontraremos una pantalla de inicio similar
a la que se presenta a continuación:
El lector recordará que usaremos Python 3, con distribución Anaconda, trabajando sobre Jupyter Notebook. En
este punto, si se han seguido las instrucciones dadas, ya contamos con la distribución Anaconda para Python 3
en nuestro computador; lo único que hace falta es tener acceso a Jupyter Notebook, lo que es, a decir verdad,
el paso más simple de todos: tan solo consiste en hacer clic.
El lector puede observar que Jupyter Notebook es una de las aplicaciones presentes en la pantalla de inicio
de Anaconda Navigator; para ingresar a ésta, tan solo debe pulsar en Launch:
6
Una vez pulsado, seremos automáticamente redirigidos a localhost. El proceso de redirección puede demorar
un poco, por lo cual, se recomienda al usuario ser paciente; una vez culminado este proceso, estaremos en
Jupyter Notebook, donde finalmente podremos iniciar con nuestro análisis econométrico y entrar en contacto
efectivo con Python.
La descripción dada hasta ahora acerca del proceso a seguir para acceder a Jupyter Notebook se ha pre-
sentado teniendo en consideración la instalación inmediatamente previa de Anaconda; sin embargo, una vez se
cuenta con ésta, no es necesario ingresar a Jupyter Notebook a través de Anaconda Navigator. Una ruta más
breve, y con una ejecución más veloz, es el acceso desde Anaconda Prompt; basta con hacer la búsqueda y
pulsar sobre la aplicación correspondiente:
7
Una vez se ha hecho clic, se debe escribir jupyter notebook en Anaconda Prompt y pulsar la tecla Enter
(Intro), después de lo cual se abrirá Jupyter Notebook y aparecerá en Anaconda Prompt información adicional
sobre el proceso ejecutado.
Por último, el modo más sencillo (una vez se tiene instalada Anaconda) y, naturalmente, el más evidente
para acceder a Jupyter Notebook es buscando directamente la aplicación y pulsando sobre el ícono correspon-
diente:
8
Cualquiera de los 3 métodos de acceso señalados anteriormente es completamente válido y debe permitir
al usuario ingresar a Jupyter Notebook; estando en esta aplicación, se debe hacer clic en la pestaña New que se
encuentra en la sección superior derecha, y pulsar entonces sobre Python 3, como se indica a continuación:
Una vez hemos hecho clic sobre Python 3, se abrirá una nueva pestaña en el navegador, que tendrá una
apariencia como la siguiente:
9
Es en este punto donde empieza nuestro verdadero trabajo, pues es ahora cuando podremos comenzar a
escribir código Python y llevar a cabo el análisis econométrico que tenemos planeado.
El lector observará que la apariencia del notebook de Jupyter es bastante agradable, con cierta similitud
a un editor de texto tradicional o un programa ofimático común. En realidad, la sensación de seguridad que
transmite esta apariencia ’familiar’ efectivamente se ve reflejada en facilidad de manejo: vemos un entorno
no sobrecargado y la mayoría de comandos se explican por sí mismos; al usuario no completamente ajeno a
la informática le resultará verdaderamente sencillo el uso del notebook de Jupyter. Tal vez el lector ha notado
que el logo de Python se encuentra en la parte superior derecha del notebook, y que cerca de éste está escrito
Python 3: esto nos indica que el código escrito en el notebook corresponde a código Python. Asimismo,
observamos que, en la parte superior del notebook, al lado del logo de Jupyter Notebook, se encuentra la
palabra Untitled; éste es el nombre que ha sido asignado por defecto a nuestro notebook, el cual, si lo
deseamos, podemos cambiar por el de nuestra preferencia haciendo clic sobre Untitled, o recurriendo a las
opciones Save as y Rename del menú File:
Modificaremos el nombre de nuestro notebook y lo llamaremos, con fines ilustrativos, “Mi primer note-
book”:
Ahora, procederemos a escribir nuestra primera línea de código; para esto, emplearemos la función prede-
finida print(...). El lector puede observar que en el notebook hay una celda sin ningún contenido; es aquí
donde debemos escribir nuestro código. Para este ejemplo específico, el argumento de la función print(...)
será “Esta es la primera línea de código de mi notebook escrito en Python”:
A esta altura, el usuario ha escrito su primera línea de código Python en el notebook, pero, ¿este código ha
tenido algún efecto? ¿se ha ejecutado alguna acción? La respuesta es simplemente no. Nuestra línea de código
10
no ha generado ningún resultado ya que tan solo la hemos escrito; aún no la hemos ejecutado. Para ejecutar
el código, debemos pulsar el botón Run, ubicado en la parte superior, debajo de la ficha Cell, o, de forma
alternativa y más rápida, emplear la combinación de teclas Ctrl + Enter.
Al ejecutar el código, obtendremos lo siguiente:
Podemos observar que debajo de nuestra celda ha sido impreso el mensaje contenido en el código y que
en los corchetes que se encuentran al lado de la celda de código, los cuales estaban vacíos, ha aparecido el
número uno. Este número indica que es la primera ejecución de código que hemos realizado en el notebook; si
ingresamos de nuevo a la celda y ejecutamos otra vez el código (con Run o la combinación Ctrl + Enter),
observaremos que ya no está el número uno sino el número dos (se invita al lector a comprobarlo por su propia
cuenta).
Ya hemos creado un notebook de Jupyter, le hemos asignado un nombre y hemos escrito y ejecutado algo
de código dentro de éste. Estas son algunas de las instrucciones básicas que el usuario, sin excusa alguna, debe
conocer; por último, guardaremos los cambios realizados en el notebook y procederemos a detenerlo y cerrarlo.
Conociendo este conjunto de operaciones, el usuario tendrá la preparación fundamental requerida para abordar
exitosamente el verdadero contenido en el que se especializa esta obra.
El proceso empleado para guardar los cambios efectuados en el notebook es muy sencillo y prácticamente
idéntico al utilizado con cualquier programa ofimático común, por lo que no representará dificultad alguna a
cualquier usuario con conocimiento básico de informática. Basta con usar la combinación de teclas Ctrl + S,
o hacer clic sobre la ficha de Save and Checkpoint (identificada con la forma de un disquete) de la cinta de
opciones.
Aunque el notebook puede cerrarse cerrando la pestaña del navegador en la que se encuentra, tal proce-
dimiento resulta no recomendable. Para cerrar correctamente el notebook, se debe pulsar sobre Close and
Halt del menú File, como se indica a continuación:
Habiendo realizado tal proceso, la pestaña del navegador se cerrará automáticamente y seremos redirigidos
a la pantalla de inicio de Jupyter Notebook, donde debemos encontrar el notebook que hemos apenas creado,
guardado y cerrado:
Podemos observar que nuestro notebook “Mi primer notebook”tiene un tamaño de 814 B y ha sido guar-
dado con extensión .ipynb; esto se debe a que los notebooks de Jupyter son archivos que se almacenan auto-
máticamente con dicha extensión. Para acceder de nuevo al notebook basta con hacer clic sobre su nombre;
de este modo, se abrirá una nueva pestaña en el navegador en la cual se encontrará el notebook.
11
Es importante resaltar que cada vez que se abre el notebook se está “iniciando una nueva sesión”, por lo que
ninguna línea de código ha sido ejecutada y, aunque es posible visualizar el output de las líneas de código, el
output que se visualiza no es el resultado de ejecuciones de la “sesión” actual. Esto es importante ya que, si se
pretende ejecutar una línea de código particular que requiere la ejecución previa de otra línea específica y esta
última no se ha ejecutado, se obtendrá un error.
Ya se ha señalado el proceso para ingresar a Anaconda Navigator y a Jupyter Notebook (mencionando 3
alternativas), ahora es momento de abordar el proceso contrario: cómo salir. Aunque es posible salir de Jupyter
Notebook con tan solo cerrar el navegador en el que se encuentra, tal proceder resulta inadecuado. El método
correcto para salir de Jupyter Notebook es cerrar (apropiadamente) los notebooks activos y, en seguida, cerrar
Jupyter Notebook haciendo clic sobre la ficha Quit que se encuentra en la sección superior derecha:
Ahora sí, podemos cerrar con entera tranquilidad el navegador. En cuanto a Anaconda Navigator, después
de haber realizado el proceso apenas descrito, si está abierta, basta con cerrarla como si fuera un programa
ofimático cualquiera, con tan solo pulsar la ficha con la x en el extremo superior derecho.
Para volver a acceder a Jupyter Notebook tan solo se requiere repetir las instrucciones necesarias que se
han indicado a lo largo de este capítulo, las cuales, en realidad, no son más que acceder a Anaconda Navigator y,
desde esta aplicación, ingresar a Jupyter Notebook (o simplemente acceder desde Anaconda Prompt o Jupyter
Notebook directamente) siendo redirigido a localhost.
Con este capítulo se cierra el tratamiento de los programas a emplear; ya es hora de entrar en la verdadera
materia que aborda esta obra, por lo que, a partir del siguiente capítulo, se desarrollará el propósito de esta
guía práctica: análisis econométrico con Python.
12
Capítulo 2
Resulta imposible llevar a cabo un análisis sin tener datos que analizar; esto no es más que un absurdo. Así,
si queremos realizar un ejercicio de análisis econométrico, debemos contar con la información requerida para
poder llevarlo a cabo; el lector no debe preocuparse por este comentario, en cuanto, además de la guía práctica
que se le ofrece, también se le brinda la información con la cual se construye los modelos y se estructura el
análisis.
Para cada uno de los ejercicios que se desarrollarán en esta obra se indicará la fuente de la cual se ha ob-
tenido la información empleada en el análisis econométrico correspondiente. Asimismo, es posible obtener los
datasets utilizados, en los mismos formatos que se emplean en esta obra, contactando al autor de la misma
(basta con enviar un correo a la dirección fatrianaa@[Link]), quien le dará acceso a la información reque-
rida, o, alternativamente, acceder al sitio web de la Unidad de Informática y Comunicaciones de la Facultad de
Ciencias Económicas de la Universidad Nacional de Colombia.
Esta obra, dado su carácter práctico y la utilidad que puede llegar a significar para los estudiantes de pre-
grado en Economía, se apoya, principalmente, en Econometría de Gujarati y Porter (2010) e Introducción a la
econometría. Un enfoque moderno de Wooldridge (2010), textos muy sencillos ampliamente empleados como
guía en los cursos básicos de Econometría, a los que en este trabajo se recurre, en múltiples ocasiones, pa-
ra desarrollar los ejercicios ilustrativos planteados y con los que se pretende que el lector pueda contrastar
resultados, de modo que evidencie que los procedimientos realizados son correctos y llegan a las soluciones
esperadas.
El dataset que se utilizará para la regresión lineal simple por Mínimos Cuadrados Ordinarios que se tratará
en el próximo capítulo es el empleado en el Ejemplo 7.1, Mortalidad infantil en relación con el PIB per cápita y
la tasa de alfabetización de las mujeres, de Gujarati y Porter (2010). Se invita al lector a consultar el ejemplo, de
modo que pueda verificar los resultados obtenidos.
La fuente de la información usada por Gujarati y Porter (2010) es “Chandan Mukherjee, Howard White y
Marc Whyte, Econometrics and Data Analysis for Developing Countries, Routledge, Londres, 1998, p. 456”. Los
datos usados para el ejercicio en Python que se lleva a cabo en esta obra se obtuvieron del fichero de datos de
Gujarati en gretl.
En este punto empieza nuestro trabajo con Python. Lo primero que debe hacer el usuario es acceder a
Jupyter Notebook (a través de Anaconda Navigator o por el medio de su preferencia) y crear un nuevo notebook.
En este notebook es donde se escribirá el código requerido y se realizará el análisis econométrico planteado.
13
In [1]: # Importamos las librerías requeridas:
import numpy as np
import pandas as pd
import [Link] as sm
import [Link] as plt
import seaborn as sns
%matplotlib inline
[Link]("seaborn-white")
¿De qué tratan estas líneas de código? El lector notará que lo primero que se ha escrito, a modo de co-
mentario (los comentarios van precedidos de #), es “Importamos las librerías requeridas”; esto es lo que, efec-
tivamente, se consigue al ejecutar esta celda de códigos. Básicamente, lo que logramos es preparar nuestro
entorno de trabajo al obtener acceso al conjunto de herramientas que requeriremos para nuestro análisis; tales
herramientas hacen parte de librerías.
Una librería simplemente es un conjunto de funciones y métodos diseñados para realizar tareas específi-
cas, agrupados en un único espacio (la librería), y que se agrega a Python “básico” con el objetivo de alcanzar
resultados que con sus funciones integradas (nativas) no pueden lograrse fácilmente. Al lector familiarizado
con el lenguaje de programación R basta con decir que una librería Python es, esencialmente, equivalente a un
paquete de R.
Las librerías que importamos son NumPy, pandas, Statsmodels, Matplotlib y Seaborn; cada una de estas tiene
un rol particular a desempeñar en nuestro trabajo, el cual se presentará, junto a breves descripciones, en un
momento. Algo que debe notar el lector es que la importación de las librerías no consiste tan solo en “importar
las librerías”, sin especificar algún tipo de instrucción adicional: la librería NumPy no se importa simplemente
como import numpy sino que se importa como import numpy as np. ¿A qué se debe esto?, ¿qué efecto
genera?
El proceso de importación con as, que es una palabra reservada de Python, se lleva a cabo con el propósito
de asignar un alias a las librerías importadas. El lector puede estar preguntándose: si la librería ya tiene un
nombre, ¿para qué se le quiere asignar un alias? La respuesta al interrogante es muy sencilla: tan solo por
practicidad. Al importar la librería NumPy con el alias de np, tal librería queda identificada con este nombre
(np); por lo tanto, cada vez que se necesite emplear una función o método de NumPy basta con hacer referencia
a np y no a numpy.
Tal vez la utilidad de as resulte más evidente en el caso de la librería Matplotlib: en vez de hacer referencia
a [Link] en cada ocasión que se requiera una función o método proveniente de éste, basta con
hacer referencia a plt. Así, en vez de usar 17 caracteres, tan solo se requiere emplear 3; aunque al lector esto le
parezca trivial en este momento, con el tiempo verá que la importación con as resulta tremendamente práctica
y facilita enormemente el trabajo.
Las librerías importadas son NumPy, pandas, Statsmodels, Matplotlib y Seaborn. ¿Por qué se importa este grupo
específico de librerías y no otras? Para comprenderlo, se ofrece una breve descripción de cada una de estas, de
modo que el lector tenga claridad sobre el propósito de las mismas.
NumPy es una librería Python para computación científica que provee un conjunto de rutinas para la
ejecución de procedimientos sobre objetos, como vectores y matrices, los cuales incluyen, entre otros,
operaciones matemáticas y lógicas, álgebra lineal básica, simulación aleatoria y estadística fundamental
(The SciPy community, 2019).
pandas es una librería de código abierto que ofrece herramientas de alto desempeño y fáciles de usar
para el manejo y análisis de datos en Python. Se trata de una librería que permite llevar a cabo un proceso
óptimo de análisis de información directamente en Python, sin tener que recurrir a un lenguaje con mayor
especificidad, como R. pandas es empleada tanto en ámbitos académicos como comerciales, en campos
14
diversos que incluyen las finanzas, la neurociencia, la estadística y la economía (Pandas, s.f). El creador
y Benevolent Dictator for Life de pandas es el matemático y desarrollador estadounidense Wes McKinney.
Statsmodels es una librería para análisis econométrico y estadístico en Python que provee funciones para
la estimación de múltiples modelos estadísticos, así como la realización de pruebas y análisis de infor-
mación estadística. Los resultados generados por las funciones de la librería son contrastados con los
de paquetes estadísticos reconocidos, de modo que se pueda garantizar su exactitud. La librería tiene
su origen en el módulo models de [Link] y en 2009, tras un proceso de corrección, prueba y
mejora, se lanzó Statsmodels de forma independiente (Perktold, Seabold, y Taylor, 2018).
Matplotlib es una librería de gráficos 2D usada en Python para la generación de imágenes de alta cali-
dad. Con Matplotlib es posible generar histogramas, gráficos de barras, diagramas de dispersión, entre
otros, mediante el uso de tan solo unas cuantas líneas de código. El módulo pyplot provee una interfaz
similar a MATLAB, útil para gráficos sencillos, que permite un control completo de las distintas propieda-
des. Matplotlib es una contribución de John Hunter, un neurobiólogo estadounidense, y de sus múltiples
colaboradores (The Matplotlib development team, 2018).
Seaborn es una librería para la elaboración de gráficos estadísticos en Python. Está basada en Matplotlib
y sus funciones de construcción de gráficos se orientan a datasets, por lo que está estrechamente inte-
grada con las estructuras de datos de pandas. El propósito de Seaborn es constituir la visualización en un
elemento clave en la exploración y comprensión de información (Waskom, 2018).
El lector, teniendo en consideración las breves descripciones apenas presentadas, ya debería tener una idea
del papel que desempeña cada una de las librerías mencionadas en el análisis que se pretende llevar a cabo; sin
embargo, si éste no es el caso y la información previamente expuesta le resulta un poco confusa, en palabras
sencillas la utilidad de cada una de las librerías es la siguiente:
Una vez ejecutadas las líneas de código de la primera celda, se habrán importado las librerías requeridas y
tendremos a disposición el conjunto de herramientas de las cuales haremos uso para llevar a cabo el análisis
econométrico. Ahora, procederemos a cargar la información con la que trabajaremos.
15
In [2]: data = pd.read_csv("C:/Users/FCE/Documents/Econometrics/[Link]",
sep = " ", delimiter="\t")
El lector debe notar que no se ha escrito read_csv(...) sino pd.read_csv(...). Esto se debe a que
read_csv(...) no es una función “nativa” de Python sino que pertenece a pandas, por lo cual es necesa-
rio especificar de dónde proviene. Dado que la importación se realizó con import pandas as pd, el código
pd.read_csv(...) informa que la función read_csv(...) pertenece a pandas. Si la importación no se hu-
biera llevado a cabo con el uso de un alias, se debería haber recurrido a pandas.read_csv(...) para conse-
guir el resultado.
Nota: El usuario debe tener en cuenta que el valor que asignará al parámetro filepath_or_buffer no será el
mismo que el del código del ejemplo anterior, pues éste depende de la ubicación del archivo GujaratiPor-
[Link] (o el archivo en el que haya guardado el dataset) en su computador.
Una vez ejecutada la línea de código de esta celda, el dataset debería haber sido importado correctamente.
Para comprobarlo, se recurre al método .head(), aplicado al objeto que contiene el dataset (en este caso se
le ha asignado el nombre data):
In [3]: [Link]()
Podemos observar que el proceso de importación ha sido exitoso y la información del dataset ha sido alma-
cenada correctamente en un DataFrame de pandas. Ahora, podemos continuar tranquilamente con el desarrollo
de nuestro análisis, en cuanto hemos concluido satisfactoriamente el primer paso.
Estadística descriptiva
Antes de empezar a examinar los valores de las variables del dataset e identificar las posibles relaciones que
hay entre estas, es recomendable hacer un reconocimiento del propio dataset, es decir, obtener información
sobre sus dimensiones, el tipo de datos que contiene, etc. Para conocer las dimensiones del dataset puede
recurrirse al atributo .shape:
In [4]: [Link]
Out[4]: (64, 4)
El resultado de la aplicación del atributo .shape es una tupla Python (...) en la que se informa el número de
filas y de columnas (en ese orden) que tiene el DataFrame; así, nuestro dataset contiene 64 filas y 4 columnas.
Ahora, sabemos que son 4 columnas, pero ¿cómo se llaman? ¿cuál es el nombre de cada una de estas 4 co-
lumnas? Cuando se ha aplicado el método .head() se han visualizado las primeras observaciones del dataset,
incluyendo el encabezado en el que se encuentran los nombres de las columnas; sin embargo, si específicamente
se quiere saber el nombre de las columnas, puede emplearse el atributo .columns:
16
In [5]: [Link]
Con toda la información obtenida se sabe que el dataset con el que se trabajará tiene 64 filas y 4 columnas,
y que los nombres de las columnas son ’CM’, ’FLR’, ’PGNP’ y ’TFR’. Sin embargo, aún se desconoce qué tipo de
datos contienen dichas columnas, ¿se trata de palabras? ¿de números? ¿solo números enteros? Para descubrirlo,
puede emplearse el atributo .dtypes:
In [6]: [Link]
Out[6]: CM int64
FLR int64
PGNP int64
TFR float64
dtype: object
La información generada por el atributo .dtypes indica que tres variables (columnas) contienen datos de
tipo int64 (número entero) y una de tipo float64 (número decimal).
Finalmente, una manera de obtener la información que se halla con la aplicación de los tres atributos apenas
descritos es el método .info(), el cual permite visualizarlos en un único espacio:
In [7]: [Link]()
<class '[Link]'>
RangeIndex: 64 entries, 0 to 63
Data columns (total 4 columns):
CM 64 non-null int64
FLR 64 non-null int64
PGNP 64 non-null int64
TFR 64 non-null float64
dtypes: float64(1), int64(3)
memory usage: 2.1 KB
El resultado de la aplicación de este método es tan solo un resumen de la información que hemos obtenido
previamente. Adicionalmente, se presenta el dato sobre el uso de memoria asociado, que en este caso es de 2.1
Kb; el tipo de objeto al que corresponde el dataset, el cual es un DataFrame de pandas; y el rango del índice,
el cual corresponde al intervalo [0, 63] (es muy importante tener en cuenta que la indización en Python inicia
desde 0 y no desde 1).
Para obtener estadística descriptiva de las variables numéricas del dataset se puede recurrir al método
.describe():
In [8]: [Link]()
17
min 12.000000 9.000000 120.000000 1.690000
25 % 82.000000 29.000000 300.000000 4.607500
50 % 138.500000 48.000000 620.000000 6.040000
75 % 192.500000 77.250000 1317.500000 6.615000
max 312.000000 95.000000 19830.000000 8.490000
Tal y como observamos, el método .describe() nos permite obtener información estadística básica sobre
las variables numéricas, lo que incluye, por ejemplo, el promedio, los valores máximo y mínimo y la desviación
estándar. Así, se puede evidenciar, por ejemplo, que el valor mínimo de la variable ’CM’ es 12, el valor máximo
de la variable ’PGNP’ es 19,830 y el promedio de la variable ’TFR’ corresponde a 5.55.
Al aplicar el método .describe() sobre el DataFrame se obtiene, por defecto, información estadística
básica de todas las variables numéricas presentes en éste. Si solo se está interesado en la información de una
variable en particular, esta puede ser seleccionada por medio del uso de los corchetes [ ], indicando su nombre
(entre comillas) dentro de estos.
Así, si, por ejemplo, solo nos interesa la variable ’CM’, el código a emplear será el siguiente:
In [9]: data["CM"].describe()
Si, por el contrario, lo que nos interesa es un grupo de variables y no solo una, podemos seleccionarlo
haciendo uso de los corchetes [ ] e incluyendo los nombres de las variables de interés (entre comillas, separados
por comas) dentro de una lista Python […]:
A esta altura se sabe que el dataset contiene información sobre 4 variables numéricas (3 de valores enteros
y 1 de valores no enteros), contando con datos de 641 individuos (entidades) para cada una de estas. Ya se ha
1 Eneste caso, la información de cada individuo (entidad) corresponde a una única fila, por lo que, como el dataset tiene 64 filas,
sabemos que hay información sobre 64 individuos (entidades). Sin embargo, en otros casos, como los paneles de datos, la información
del mismo individuo (entidad) no corresponde a una única fila, por lo que hay que ser cuidadosos en tales situaciones.
18
realizado un muy breve reconocimiento de los datos, por lo que se tiene un mayor nivel de familiarización con
los mismos y es posible identificar potenciales relaciones a examinar.
Hasta el momento, se ha llevado a cabo un proceso de reconocimiento de los datos por medio del cual se
ha obtenido información puramente numérica: el número de filas y columnas, la media y la desviación estándar
de cada variable, etc. Esto nos da cierta idea sobre los datos, sin embargo, con frecuencia, “una imagen vale más
que mil palabras”2 ; es posible que un gráfico tenga un inmenso poder comunicativo y contribuya enormemente
a obtener una mejor comprensión acerca de los datos con los que se trabajará.
La función pairplot(...) de Seaborn resulta particularmente útil para el contexto en el que nos encontra-
mos; esta función crea una cuadrilla en la que podremos visualizar la relación que existe entre pares de variables
(por medio de diagramas de dispersión) y la distribución de cada una de estas (por medio de histogramas):
19
En la diagonal se observa la distribución de cada una de las variables, mientras que en los diagramas de
dispersión se encuentra la relación entre pares de estas; así, se puede observar que, por ejemplo, el ingreso
se concentra en valores inferiores a 5000 y existe una relación negativa entre la tasa de alfabetización de las
mujeres y la mortalidad infantil.
En estos momentos ya se cuenta con información básica sobre las variables del dataset y hemos adelantado
una breve inspección visual de las relaciones que se presentan entre estas; ahora, se puede hacer uso de este
conocimiento para obtener expresiones más concretas de dichas relaciones, examinándolas desde un punto de
vista un tanto más técnico: en el siguiente capítulo se abordará el tema de la regresión lineal por el método de
Mínimos Cuadrados Ordinarios.
20
Capítulo 3
La regresión lineal es el tema introductorio que generalmente se aborda en los cursos de Econometría de
pregrado en Economía; de acuerdo a Greene (2003), “el modelo de regresión lineal es la herramienta más útil del
kit de herramientas de los econometristas” (p. 7). Este tema suele ser el primer contacto que tiene el estudiante
con materia de contenido propiamente econométrico y tiende a ser abordada con cierto detalle, de modo que
éste tenga claridad suficiente sobre la misma. Es común que en los cursos se presente la regresión lineal, junto
a los fundamentos teóricos aplicables a la misma, con el apoyo de un texto guía de un autor reconocido y
que se haga énfasis apreciable en los supuestos asociados, de modo que el estudiante adquiera un nivel de
conocimiento adecuado sobre ésta.
En esta obra, considerando su carácter práctico y teniendo en cuenta la advertencia realizada en el prefacio
de que los fundamentos teóricos se abordarán con un nivel de profundidad mínimo, se estará limitado a exponer
la idea general de la regresión lineal en términos muy superficiales; esto, con el propósito de enfocarse en el
objetivo planteado y de no hacer sentir incómodo al lector proveniente de otras ramas de estudio o con una
formación econométrica no tan sólida, al no ahondar en las bases técnicas y las formalidades matemáticas sobre
las que se cimienta la regresión lineal y sus métodos de estimación.
El análisis de regresión trata del estudio de la dependencia de una variable (variable dependiente)
respecto de una o más variables (variables explicativas) con el objetivo de estimar o predecir la
media o valor promedio poblacional de la primera en términos de los valores conocidos o fijos (en
muestras repetidas) de las segundas. (Gujarati y Porter, 2010, p.15)
En términos simples, el análisis de regresión estudia la relación que existe entre la variable dependiente y
la(s) variable(s) explicativa(s). Para el caso de la regresión lineal, tal relación puede ser expresada como una
función lineal en los parámetros, pero ¿de qué parámetros se está hablando?
Para resolver el anterior interrogante, lo primero que se hará será considerar una regresión bivariada, en la
cual la variable dependiente (que se identifica con y) puede ser expresada como una función1 (con linealidad
en los parámetros) de la variable explicativa (que se identifica con x). Concretamente:
yi = β 0 + β 1 xi + ui
Esta expresión corresponde a la Función de Regresión Poblacional (FRP) (Gujarati y Porter, 2010) y los
parámetros a los que se ha hecho referencia son simplemente los β; como podemos ver, estos β no se ven
afectados por algún tipo de transformación que modifique su relación lineal con la variable dependiente: a esto
1 No se trata exactamente de una función, pues, como lo señalan Fahrmeier, Kneib, y Lang (2007), una característica fundamental
de los análisis de regresión es que la relación entre la variable dependiente y las variables explicativas no corresponde a una función en
sentido estricto, dado que se ve afectada por perturbaciones aleatorias (p. 19).
21
es que se hace referencia con linealidad en los parámetros, la cual, no necesariamente debe aplicar a las variables
explicativas (Greene, 2003; Gujarati y Porter, 2010). En cuanto a u, ésta corresponde al término de error, el cual
recoge el efecto de las variables no incluidas explícitamente en el modelo y que afectan y.
Ahora, en la práctica se desconoce la Función de Regresión Poblacional, por lo que se recurre a la Función
de Regresión Muestral (FRM) como una aproximación a esta. Al lector interesado en la cuestión y que desee
conocer una explicación al respecto, se sugiere consultar el Capítulo 2, Análisis de regresión con dos variables:
algunas ideas básicas, de “Gujarati, D.N. y Porter, D.C. (2010). Econometría”.
La Función de Regresión Muestral es:
yi = β̂ 0 + β̂ 1 xi + ûi
Esta expresión es muy parecida a la de la Función de Regresión Poblacional, pero presenta una diferencia
fundamental: los términos de la Función de Regresión Muestral tienen correspondencia directa con los términos
incluidos en la expresión de la Función de Regresión Poblacional, pero no son exactamente los mismos, en
cuanto se trata de sus estimadores. Así, el estimador del parámetro k se identifica como k̂.
Tal y como lo señalan Gujarati y Porter (2010), “el objetivo principal del análisis de regresión es estimar la
Función de Regresión Poblacional con base en la Función de Regresión Muestral” (p. 44). En este sentido, lo
que se busca es hallar los valores de los β̂ que sean más cercanos a los verdaderos valores de los β.
Para realizar la estimación, es frecuente el uso de dos métodos: Mínimos Cuadrados Ordinarios y Máxima
Verosimilitud. El método de Mínimos Cuadrados Ordinarios, o MCO, es el más común y tiende a ser abordado
como tema fundamental en los cursos de econometría a nivel introductorio.
Resaltando, una vez más, que el nivel de profundidad con el que se abordan las bases teóricas en esta obra
es ínfimo, tan solo se presentará (con el perdón de los especialistas, al no tratar con el detalle justo la cuestión)
la idea general del método de MCO: el objetivo es, básicamente, minimizar la suma de los residuos (termino û
en la Función de Regresión Muestral) cuadrados.
Para el caso de una regresión bivariada, lo que se busca con el método MCO es el menor valor posible de
n
∑i=1 û2i , teniendo en cuenta que:
n n
∑ û2i = ∑ (yi − β̂0 − β̂1 xi )2
i =1 i =1
Así, los estimadores β̂ 0 y β̂ 1 con los que se halla la menor suma de los residuos al cuadrado son los estima-
dores MCO, los cuales, con el cumplimiento de un conjunto de supuestos específicos, evidencian propiedades
estadísticas muy atractivas.
Al lector interesado en conocer el proceso por medio del cual se obtienen los valores de los β̂ y familiarizarse
con las bases matemáticas correspondientes, se le invita a consultar el Capítulo 3, Modelo de regresión con dos
variables: problema de estimación, de “Gujarati, D.N. y Porter, D.C. (2010). Econometría” y el Capítulo 2, El modelo
de regresión simple, de “Wooldridge, J.M. (2010). Introducción a la econometría. Un enfoque moderno”.
Ya se ha presentado, muy brevemente, una idea general de la regresión lineal y el método de MCO. Ahora, es
momento de iniciar de forma efectiva con la labor práctica; a saber, análisis econométrico con el uso de Python.
Lo primero que se hará es llevar a cabo una estimación, por Mínimos Cuadrados Ordinarios, de un modelo
de regresión lineal simple. El modelo de regresión lineal simple no es más que una regresión lineal bivariada,
es decir, en la que solo intervienen dos variables explícitas: la variable dependiente como función lineal en los
parámetros de la variable explicativa. Esto es, simplemente, un modelo al que corresponde una FRM de la forma
yi = β̂ 0 + β̂ 1 xi + ûi , la cual debe resultar familiar al lector, pues es la que se ha tratado previamente.
22
Regresión lineal simple
Para llevar a cabo la regresión lineal simple, utilizaremos el mismo dataset del capítulo previo y reque-
riremos de las mismas herramientas empleadas en éste, por lo que procederemos a importar las respectivas
librerías y los datos:
En el capítulo previo se llevó a cabo un breve reconocimiento de los datos. Ahora, se procederá a la reali-
zación de una Regresión Lineal Simple por el método de Mínimos Cuadrados Ordinarios, teniendo en conside-
ración lo siguiente:
En la regresión lineal simple que se adelantará, la variable dependiente será ’CM’ (mortalidad infantil) y
la variable explicativa será ’PGNP’ (PIB per cápita).
Se procederá a examinar visualmente la relación que existe entre las variables de interés. Para esto, se
construirá un diagrama de dispersión (scatterplot) con el uso del método .scatter(), aplicado a un Axes de
Matplotlib. El código a emplear es el siguiente:
23
fontsize = 13, fontweight = "bold",
ha = "right")
[Link](.9,-.08,
"Triana, F.\n(2019)",
fontsize = 12, ha = "right")
[Link]()
Aunque es un bloque de código algo extenso, el lector debe centrar su atención tan solo en la primera,
especialmente la tercera, y la última línea de código de la celda, pues son las que contienen la “esencia” del
gráfico; las demás, solamente hacen referencia a detalles de los cuales se puede prescindir. No se hará énfasis
en la explicación de este código, en cuanto no es el asunto principal que busca tratarse y la mayoría de líneas “se
explican por sí mismas”. Al lector que, sin embargo, esté interesado en este asunto, se recomienda el estudio
detallado de Matplotlib.
A partir del gráfico se puede observar que la relación entre las variables es aparentemente negativa. Sin
embargo, esto tan solo es una impresión que se genera a partir de la observación y puede ser errada, razón por
la cual se examinará cuantitativamente dicha relación buscando comprobar si la conjetura es correcta; para tal
propósito se llevará a cabo una regresión lineal por el método de MCO.
Para crear el modelo, se emplea la función OLS(...) de Statsmodels. Los argumentos que se incluyen co-
rresponden, respectivamente, a la variable dependiente y a la variable explicativa. El resultado será asignado a
un objeto que se denominará MiModeloSimple:
El modelo ya ha sido creado, es decir, ya se ha definido su estructura; sin embargo, aún no se ha llevado
a cabo ninguna estimación a partir de éste. Para efectuar la estimación se recurrirá al método .fit() y se
almacenarán los resultados generados en un objeto al que se denominará ResultadosSimple:
24
In [5]: ResultadosSimple = [Link]()
Al ejecutar esta línea de código no se obtiene algún resultado visible. Esto se debe a que lo que se ha
conseguido con la ejecución es guardar los resultados generados en un objeto, sin indicar la realización de
alguna acción en particular con dicho objeto. Ahora, si se desea visualizar los resultados del modelo, puede
simplemente usarse la función integrada print(...), empleando como argumento la aplicación del método
.summary() sobre el objeto que almacena los resultados de dicho modelo. Así:
In [6]: print([Link]())
OLS Regression Results
==============================================================================
Dep. Variable: CM R-squared: 0.056
Model: OLS Adj. R-squared: 0.041
Method: Least Squares F-statistic: 3.710
Date: xxx, xx xxx xxxx Prob (F-statistic): 0.0586
Time: xx:xx:xx Log-Likelihood: -413.92
No. Observations: 64 AIC: 829.8
Df Residuals: 63 BIC: 832.0
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
PGNP 0.0124 0.006 1.926 0.059 -0.000 0.025
==============================================================================
Omnibus: 7.668 Durbin-Watson: 0.755
Prob(Omnibus): 0.022 Jarque-Bera (JB): 7.941
Skew: -0.561 Prob(JB): 0.0189
Kurtosis: 4.312 Cond. No. 1.00
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
Lo que se obtiene con la ejecución del código es un resumen de la instancia de resultados en el que se
presenta la información más relevante acerca de estos: en la parte superior se incluye el método de estimación
empleado, el número de observaciones, los grados de libertad (de los residuos y del modelo) y el R2 (estándar
y ajustado), entre otros datos relevantes. En la sección del medio se encuentran los coeficientes, junto a sus
respectivos errores estándar, estadísticos t e intervalos de confianza. En la última sección están las advertencias,
las cuales resultan, frecuentemente, muy útiles y pueden ayudar a identificar potenciales problemas.
Se puede observar que el coeficiente de la variable explicativa es positivo, lo que contradice la conclusión
visual de que la relación, aparentemente, era negativa; algo que resulta, además, poco lógico, pues se esperaría
que los países con ingreso per cápita más alto tuvieran una tasa de mortalidad infantil menor. Ante esta situa-
ción extraña vale la pena preguntarse sobre su causa: ¿cómo puede explicarse tal inconsistencia? La respuesta
es extremadamente simple: el lector debe notar que para este modelo solo se ha calculado el coeficiente co-
rrespondiente a la variable ’PGNP’ y no se ha considerado la existencia del parámetro de intercepto, es decir, se
ha llevado a cabo una regresión a través del origen.
25
En vez de considerar yi = β 0 + β 1 xi + ui , se ha asumido que la FRP es de la forma yi = β 1 xi + ui .
La función OLS(...) de Statsmodels no asume por defecto que el modelo a estimar incluye una constante
como variable explicativa y, por tanto, no realiza la estimación del coeficiente correspondiente al parámetro de
intercepto. Para conseguir la estimación de dicho parámetro es necesario especificar que una constante debe
añadirse al conjunto de variables explicativas del modelo; esto se logra, por ejemplo, empleando la función
add_constant(...) de Statsmodels, tomando como argumento el objeto correspondiente a las variables ex-
plicativas incluidas. Así:
Para realizar la estimación y visualizar los resultados correspondientes, basta con repetir la misma estruc-
tura de los códigos empleados previamente (en el caso del modelo sin término del intercepto). Así:
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
[2] The condition number is large, 3.43e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Ahora, se puede observar que el coeficiente correspondiente a la variable explicativa es negativo, tal y como
se había presumido en el examen visual, y que ésta es estadísticamente significativa. Lo que se hará enseguida
es graficar los valores originales y los valores estimados por el modelo, de modo que puedan contrastarse y
26
sea posible un examen visual de la bondad de ajuste de éste. Para esto, nuevamente se hará uso del método
.scatter(...) aplicado sobre un Axes de Matplotlib.
In [9]: fig, ax = [Link](figsize = (10,5))
[Link]("Regresión lineal simple (MCO)", fontsize = 18,
fontweight = "bold")
[Link](data["PGNP"], data["CM"], s = 50,
label = "Valores originales",
color = "darkblue")
[Link](data["PGNP"], [Link](), s = 50,
label = "Valores estimados",
color = "red")
ax.set_xlabel("PIB per cápita", fontsize = 15)
ax.set_ylabel("Mortalidad infantil\n(por cada 1000 nacidos vivos)",
fontsize = 15)
[Link](frameon = True, loc = "lower left")
plt.subplots_adjust(top=0.9)
plt.tick_params(labelsize = 15)
[Link](.9,-.02,
"Elaboración:",
fontsize = 13, fontweight = "bold",
ha = "right")
[Link](.9,-.08,
"Triana, F.\n(2019)",
fontsize = 12, ha = "right")
[Link]()
27
Se observa que la bondad de ajuste del modelo no es la mejor. ¿A qué se debe este resultado desalentador?
Si el lector examina el gráfico anterior podrá notar que la parte correspondiente a los datos originales presenta
un muy ligero parecido con el gráfico de una función de la forma y = f ( x ) = xa , con a > 0, x > 0 y
representación del valor de x en el eje de las abscisas y del valor de y en el eje de las ordenadas. Este parecido
sugiere que la forma funcional yi = β 0 + β 1 xi + ui tal vez no sea la más adecuada para modelar una relación
con estos datos, en cuanto y puede asemejarse más a una función de x de la forma f ( x ) = xa ( a 6= 0) que a
una de la forma f ( x ) = ax ( a 6= 0). Teniendo en cuenta esto, tal vez la forma yi = β 0 + β 1 1x + ui genere
mejores resultados.
Nota: Un modelo con forma funcional yi = β 0 + β 1 1x + ui es denominado modelo recíproco. (Se in-
vita al lector a consultar el Capítulo 6, Extensiones del modelo de regresión lineal con dos variables, de Gujarati
y Porter (2010), para obtener mayor información al respecto).
Para trabajar con un modelo recíproco se debe tener en cuenta que la estimación no se hará a partir de los
valores de x sino de los valores de 1x . Evidentemente, esto implica que el modelo es no lineal en la variable
x, pero, ¿esto significa que lo que se llevará a cabo es una regresión no lineal? La respuesta es no. Aunque el
modelo ya no es lineal en la variable x, el lector debe recordar que la linealidad considerada hace referencia a
los parámetros y no a las variables, por lo cual, a pesar de la no linealidad de la variable, sigue tratándose de un
modelo de regresión lineal, en cuanto la linealidad en los parámetros no se ha visto afectada.
Se ha creado el modelo tomando ’PGNP’ como variable explicativa, ahora, ésta no será incluida, sino que se
recurrirá a su recíproco: ’1/PGNP’. Los valores de ’PGNP’ están incluidos en el dataset que se ha importado, sin
embargo, los valores de ’1/PGNP’ son desconocidos. ¿Cómo se puede emplear una variable de cuyos valores no
se tiene información?
Aunque se desconocen los valores de ’1/PGNP’, aún se tiene la posibilidad de obtenerlos en tanto los valores
de ’PGNP’ son conocidos y tan solo se debe realizar una operación matemática sencilla. El cálculo de estos
valores, a pesar de ser simple, puede resultar tedioso dada la cantidad de los mismos y puede llevar a cometer
errores si se ejecuta de forma “manual”; por fortuna, pandas resulta muy útil en este contexto, pues evita la
realización de tal labor “manual” para el cálculo de cada uno de los valores, al permitir aplicar operaciones
matemáticas básicas de Python que permiten obtener el resultado fácilmente.
Se creará la variable correspondiente al recíproco de la variable original, ya que ésta no hace parte del
dataset. Para crear dicha variable, simplemente se realiza una asignación a una Series de pandas, señalando
el nombre del DataFrame correspondiente y, a continuación, especificando el nombre que quiere darse a la
variable nueva (escrito entre comillas) dentro de corchetes [ ]. Así:
In [10]: data["1/PGNP"] = 1/data["PGNP"]
Se verifica la correcta ejecución del proceso, recurriendo al método .head():
In [11]: [Link]()
Out[11]: CM FLR PGNP TFR 1/PGNP
0 128 37 1870 6.66 0.000535
1 204 22 130 6.15 0.007692
2 202 16 310 7.00 0.003226
3 197 65 570 6.25 0.001754
4 96 76 2050 3.81 0.000488
Como el lector puede observar, la última columna lleva el nombre que se ha asignado a la nueva variable
(“1/PGNP”) cuyos valores son los resultados de tomar, fila por fila, el número uno y dividirlo por el valor corres-
pondiente de la columna ’PGNP’. Para comprobar la exactitud de la operación realizada, siéntase en la libertad
de tomar la información de una fila cualquiera y ejecutar el cálculo correspondiente.
28
Ahora, estimaremos el modelo yi = β 0 + β 1 x1i + ui , teniendo en cuenta que 1x corresponde a la variable
recién creada ’1/PGNP’, y observaremos si se genera una mejora en el ajuste respecto al obtenido con la estima-
ción de yi = β 0 + β 1 xi + ui . Para la construcción y estimación del modelo y visualización de los resultados,
emplearemos, nuevamente, la función OLS(...) de Statsmodels y los métodos .fit() y .summary():
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
Ahora, procedemos a graficar los valores originales y los valores estimados, usando el método
.scatter(...) de un Axes de Matplotlib, de modo que podamos observar si efectivamente se genera una
mejora en el ajuste:
29
label = "Valores estimados",
color = "red")
ax.set_xlabel("PIB per cápita", fontsize = 15)
ax.set_ylabel("Mortalidad infantil\n(por cada 1000 nacidos vivos)",
fontsize = 15)
[Link](frameon = True, loc = "lower right")
plt.subplots_adjust(top=0.9)
plt.tick_params(labelsize = 15)
[Link](.9,-.02,
"Elaboración:",
fontsize = 13, fontweight = "bold",
ha = "right")
[Link](.9,-.08,
"Triana, F.\n(2019)",
fontsize = 12, ha = "right")
[Link]()
Como podemos observar, se presenta una mejora clara, en cuanto el ajuste parece ser más adecuado; la me-
jora es visualmente evidente, sin embargo, también se ha presentado un incremento apreciable del coeficiente
de determinación, o R2 , el cual ha pasado de 0.166 a 0.459, como el lector puede comprobar en las tablas de
resumen de cada regresión. Así, para este caso particular, la forma funcional del modelo recíproco parece ser
más apropiada que la forma funcional tradicional.
30
Regresión lineal múltiple
Ya hemos realizado una regresión lineal usando una única variable explicativa; ahora, emplearemos varias
explicativas. En un modelo de regresión lineal múltiple se tiene que la FRP es de la forma yi = β 0 + β 1 x1i +
... + β k xki + ui , de modo que se incluyen k variables explicativas (sin contar la constante de β 0 ) y se deben
estimar k + 1 parámetros (los de las variables y el término del intercepto).
En el modelo de regresión lineal múltiple que emplearemos, la variable dependiente será ’CM’ (mortalidad
infantil) y las variables explicativas serán ’PGNP’ (PIB per cápita) y ’FLR’ (tasa de alfabetización de las muje-
res). Lo primero que haremos es examinar gráficamente la relación que existe entre las variables explicativas
y la variable dependiente; esta vez, recurriremos a la función pairplot(...) de Seaborn, en cuanto esta nos
permite graficar relaciones de pares en conjuntos de datos.
La función pairplot(...) de Seaborn incluye múltiples argumentos, sin embargo, para nuestros propósi-
tos, además de (obviamente) el parámetro obligatorio data, los argumentos importantes son x_vars y y_vars; el
primero permite especificar las variables que se quiere graficar en el eje de las abscisas y el segundo las que se
quiere graficar en el eje de las ordenadas. Para nuestro caso particular, y_vars = 'CM' y x_vars =['PGNP',
'FLR'] (el lector debe notar que, dado que se trata de varias variables, deben incluirse dentro de una lista Pyt-
hon [...]).
31
A partir de nuestro gráfico podemos observar que la relación que existe entre las variables explicativas
y la variable dependiente es, aparentemente, negativa para ambas. Ahora, es momento de cuantificar tales
relaciones. Debemos asignar las variables a objetos, de modo que resulte menos tediosa la construcción del
modelo y contemos con un código más “limpio”. Para seleccionar variables de un DataFrame se deben emplear
los corchetes ’[ ]’, teniendo en consideración lo siguiente:
Cuando se trata de una sola variable, basta con escribir su nombre (entre comillas) dentro de los corche-
tes.
Cuando se trata de varias variables, se debe escribir sus nombres (entre comillas, separados por comas)
dentro de una lista Python [...] y emplear tal lista dentro de los corchetes.
La variable dependiente será asignada al objeto Y y las variables explicativas serán asignadas al objeto X,
así:
In [15]: Y = data["CM"]
X = data[["PGNP", "FLR"]]
Como ya hemos asignado las variables a los objetos correspondientes, ahora podemos emplear tales objetos
como argumentos de la función OLS(...) de Statsmodels para la construcción del modelo, recordando que sus
argumentos son, respectivamente, la variable dependiente y las variables explicativas:
Llevamos a cabo la estimación, con el método .fit(), y guardamos los resultados generados en un objeto
que se denominará Resultados:
32
Finalmente, visualizamos los resultados empleando la función integrada (nativa) print(...), usando co-
mo argumento la aplicación del método .summary():
In [18]: print([Link]())
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
Se puede evidenciar que no se ha estimado el término del intercepto, es decir, se ha efectuado una regresión
a través del origen. Esto, como se mencionó con anterioridad, se debe a que la función OLS(...) de Statsmodels
no asume por defecto que se debe incluir una constante en el conjunto de variables explicativas. Para que la
estimación se realice con un intercepto es necesario especificar que una constante debe ser empleada como re-
gresora; para tal propósito, se puede emplear la función add_constant(...) de Statsmodels, tomando como
argumento el objeto que contiene las demás variables explicativas.
La estimación y la visualización de resultados se llevan a cabo empleando los métodos que se han usado
con anterioridad (.fit() y .summary()):
33
OLS Regression Results
==============================================================================
Dep. Variable: CM R-squared: 0.708
Model: OLS Adj. R-squared: 0.698
Method: Least Squares F-statistic: 73.83
Date: xxx, xx xxx xxxx Prob (F-statistic): 5.12e-17
Time: xx:xx:xx Log-Likelihood: -328.10
No. Observations: 64 AIC: 662.2
Df Residuals: 61 BIC: 668.7
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 263.6416 11.593 22.741 0.000 240.460 286.824
PGNP -0.0056 0.002 -2.819 0.006 -0.010 -0.002
FLR -2.2316 0.210 -10.629 0.000 -2.651 -1.812
==============================================================================
Omnibus: 0.732 Durbin-Watson: 2.186
Prob(Omnibus): 0.693 Jarque-Bera (JB): 0.559
Skew: 0.228 Prob(JB): 0.756
Kurtosis: 2.949 Cond. No. 6.77e+03
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
[2] The condition number is large, 6.77e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
El lector puede observar que los resultados que hemos alcanzado han sido obtenidos con el empleo de
métodos, los cuales son, esencialmente, funciones asociadas a un objeto específico, aplicables directamente
sobre éste.
.fit() y .summary() son métodos aplicables a objetos generados por la función OLS(...), sin embargo,
estos no son los únicos métodos propios de tal tipo de objetos. Para conocer los métodos y atributos de un
objeto en Python, se emplea la función integrada (nativa) dir(...), la cual toma como argumento un objeto
de Python y devuelve sus atributos y métodos dentro de una lista Python ’[…]’.
Para conocer los métodos y atributos del objeto Resultados2, el código a emplear es el siguiente:
In [21]: dir(Resultados2)
Out[21]: ['HC0_se',
'HC1_se',
'HC2_se',
'HC3_se',
'_HCCM',
'__class__',
34
'__delattr__',
'__dict__',
'__dir__',
'__doc__',
'__eq__',
'__format__',
'__ge__',
'__getattribute__',
'__gt__',
'__hash__',
'__init__',
'__init_subclass__',
'__le__',
'__lt__',
'__module__',
'__ne__',
'__new__',
'__reduce__',
'__reduce_ex__',
'__repr__',
'__setattr__',
'__sizeof__',
'__str__',
'__subclasshook__',
'__weakref__',
'_cache',
'_data_attr',
'_get_robustcov_results',
'_is_nested',
'_wexog_singular_values',
'aic',
'bic',
'bse',
'centered_tss',
'compare_f_test',
'compare_lm_test',
'compare_lr_test',
'condition_number',
'conf_int',
'conf_int_el',
'cov_HC0',
'cov_HC1',
'cov_HC2',
'cov_HC3',
'cov_kwds',
'cov_params',
'cov_type',
'df_model',
35
'df_resid',
'diagn',
'eigenvals',
'el_test',
'ess',
'f_pvalue',
'f_test',
'fittedvalues',
'fvalue',
'get_influence',
'get_prediction',
'get_robustcov_results',
'initialize',
'k_constant',
'llf',
'load',
'model',
'mse_model',
'mse_resid',
'mse_total',
'nobs',
'normalized_cov_params',
'outlier_test',
'params',
'predict',
'pvalues',
'remove_data',
'resid',
'resid_pearson',
'rsquared',
'rsquared_adj',
'save',
'scale',
'ssr',
'summary',
'summary2',
't_test',
't_test_pairwise',
'tvalues',
'uncentered_tss',
'use_t',
'wald_test',
'wald_test_terms',
'wresid']
Cada uno de los objetos de esta lista Python [...] es un método o atributo del modelo que hemos construido.
Así, si deseamos conocer el coeficiente de determinación múltiple o R2 , el cual se identifica con el atributo
.rsquared, podemos emplear el siguiente código:
36
In [22]: print("El R2 del modelo es:", [Link])
Si estamos interesados en los coeficientes estimados, podemos recurrir al atributo .params, así:
Asimismo, los métodos permiten ejecutar acciones específicas sobre el objeto. A modo de ejemplo, el méto-
do .get_robustcov_results() permite obtener una nueva instancia de resultados asumiendo por defecto
covarianza robusta.
Warnings:
[1] Standard Errors are heteroscedasticity robust (HC1)
37
[2] The condition number is large, 6.77e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Se invita al lector a explorar los demás métodos y atributos disponibles, de modo que pueda conocer toda
la información adicional con la que puede profundizar su análisis.
En este punto se cierra el tercer capítulo de esta obra; el lector ya ha tenido un acercamiento a la regresión
lineal simple y múltiple en Python con estimación por Mínimos Cuadrados Ordinarios; asimismo, ha observado
cómo llevar a cabo un análisis exploratorio mínimo y cómo elaborar gráficos básicos. Ahora, tal y como se
mencionó en este capítulo, dados ciertos supuestos, los estimadores de MCO exhiben propiedades estadísticas
atractivas; tales supuestos corresponden al tema abordado en el siguiente capítulo.
38
Capítulo 4
Verificación de supuestos
Los estimadores de MCO presentan una serie de propiedades estadísticas atractivas dados unos supuestos;
tal afirmación se sustenta en el Teorema de Gauss-Markov, el cual, siguiendo a Gujarati y Porter (2010), puede
expresarse como: “Dados los supuestos del modelo clásico de regresión lineal, los estimadores de mínimos
cuadrados, dentro de la clase de estimadores lineales insesgados, tienen varianza mínima, es decir, son MELI”
(p. 72).
MELI (BLUE en inglés) hace referencia a “Mejores Estimadores Lineales Insesgados” y es la calificación que
adquieren los estimadores de MCO cuando se cumple con los supuestos del modelo clásico de regresión lineal,
también conocido como modelo de Gauss o modelo estándar de regresión lineal.
El modelo clásico de regresión lineal plantea, de acuerdo a Gujarati y Porter (2010), los siguientes supues-
tos:
En este capítulo se verifica el cumplimiento de algunos de estos supuestos con el uso de diversas herramien-
tas. Algo que debe tener en consideración el lector es que algunos de los supuestos anteriormente señalados
corresponden específicamente a la estructura (ecuación) del modelo (por ejemplo, el supuesto 1) mientras que
otros hacen referencia explícita a las características del término de error (por ejemplo, el supuesto 4 y el su-
puesto 5). Así, para tener mayor claridad, primero examinaremos los supuestos sobre la estructura del modelo
y posteriormente abordaremos los supuestos que se refieren al término de error, por lo que no se seguirá el
mismo orden en el que los supuestos han sido enunciados.
39
Supuestos sobre la estructura del modelo
El modelo clásico de regresión lineal asume el cumplimiento de una serie de supuestos sobre la estructura
(ecuación) del modelo. Es posible verificar el cumplimiento de algunos de estos supuestos en Python utilizando
funciones propias de algunas librerías especializadas como Statsmodels.
Estimaremos un modelo de regresión lineal múltiple y verificaremos el cumplimiento de los supuestos sobre
la estructura del mismo. El dataset que utilizaremos será, al igual que en el anterior capítulo, el del Ejemplo
7.1, Mortalidad infantil en relación con el PIB per cápita y la tasa de alfabetización de las mujeres, de Gujarati y Porter
(2010).
El lector ya debe saber que el anterior bloque de código permite cargar las herramientas necesarias para
llevar a cabo las acciones requeridas en el análisis econométrico planteado; en caso de que no tenga claridad
al respecto, se le pide consultar la sección “Preparación del entorno” del anterior capítulo, en donde se brinda
una explicación sobre tal cuestión.
El código de esta celda es exactamente igual al empleado en la preparación del entorno del capítulo anterior,
exceptuando por lo siguiente: utilizaremos, adicionalmente, el módulo stats de la librería Statsmodels, el cual
se importará con el alias de sms.
La importación de [Link] es fundamental para nuestra labor, pues es donde están conte-
nidas la mayoría de funciones que emplearemos, sin las cuales no sería posible verificar fácilmente el cumpli-
miento de los supuestos del modelo clásico de regresión lineal. El usuario debe ser cuidadoso de, en caso de
reutilizar el bloque de código de preparación del entorno del capítulo anterior, incluir correctamente la línea de
importación de [Link].
Al lector que no entienda esta línea de código, se sugiere, nuevamente, dirigirse al Capítulo 2, en donde
encontrará una explicación que seguramente aclarará sus dudas, o al Capítulo 1 de Triana y Galindo (2019).
40
Una vez ejecutada la línea de código de esta celda, el dataset debería haber sido importado correctamente.
Para comprobarlo, recurrimos al método .head(), aplicado al objeto que contiene el dataset (en este caso le
hemos asignado el nombre data):
In [3]: [Link]()
Podemos observar que el proceso de importación ha sido exitoso y la información del dataset ha sido alma-
cenada correctamente en un DataFrame de pandas. Ahora, podemos continuar tranquilamente con el desarrollo
de nuestro análisis, en cuanto hemos concluido satisfactoriamente el primer paso. Procederemos a verificar el
cumplimiento de algunos de los supuestos sobre la estructura del modelo.
In [4]: [Link]
Out[4]: (64, 4)
El resultado que obtenemos es una tupla Python (...) que informa que nuestro DataFrame es de tamaño
64 × 4, por lo que posee 64 filas y 4 columnas. Considerando que cada fila corresponde a una observación y
cada columna a una variable, el número máximo de parámetros a estimar es 63.
¿Qué sucede si el número de observaciones es inferior al número de parámetros a estimar? “Si hay menos
de k observaciones entonces X no puede ser de rango completo” (Greene, 2003, p. 14). Debido a esto, cuando
el número de parámetros a estimar excede el número de observaciones disponibles para hacer la estimación,
se obtiene que la varianza es infinita, por lo cual el método MCO no puede emplearse (James, Witten, Hastie, y
Tibshirani, 2013).
A modo de ejemplo, se realizará una regresión de prueba, con fines puramente ilustrativos, en la que se
usarán solo 3 observaciones para estimar 4 parámetros. Para esto, se seleccionarán únicamente las tres primeras
observaciones del DataFrame recurriendo a los corchetes [ ]:
41
OLS Regression Results
==============================================================================
Dep. Variable: CM R-squared: 1.000
Model: OLS Adj. R-squared: nan
Method: Least Squares F-statistic: 0.000
Date: xxx, xx xxx xxxx Prob (F-statistic): nan
Time: xx:xx:xx Log-Likelihood: 66.160
No. Observations: 3 AIC: -126.3
Df Residuals: 0 BIC: -129.0
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 3.2453 inf 0 nan nan nan
FLR 2.0501 inf 0 nan nan nan
PGNP -0.0692 inf -0 nan nan nan
TFR 26.7721 inf 0 nan nan nan
==============================================================================
Omnibus: nan Durbin-Watson: 0.854
Prob(Omnibus): nan Jarque-Bera (JB): 0.511
Skew: -0.678 Prob(JB): 0.774
Kurtosis: 1.500 Cond. No. 753.
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
[2] The input rank is higher than the number of observations.
Como el lector puede observar, los errores estándar son infinitos, y los p-values e intervalos de confianza no
están definidos, por lo que no resulta útil emplear los estimadores de MCO cuando el número de observaciones
es inferior al número de parámetros que se pretende estimar.
Ciertamente, el ejemplo que se acaba de presentar es muy sencillo; sin embargo, es posible que el dataset
a emplear tenga un gran número de variables explicativas, incluso muy superior a la cantidad de observaciones
disponibles, por lo cual el Supuesto #6, a pesar de su apariencia inocente, no debe ser tomado a la ligera.
Cuando existen muchas más variables explicativas que observaciones, se puede recurrir a diversas técnicas
y se pueden llevar a cabo procesos de selección de variables; sin embargo, tales procedimientos exceden el
alcance de esta obra y no serán abordados en la misma.
42
si el valor de una variable explicativa no cambia, resulta difícil identificar la manera como está asociada con la
variable dependiente.
Identificar el cumplimiento de este supuesto es bastante sencillo: “si la desviación estándar muestral de las
xi es cero, entonces el supuesto no se satisface; si no es así, este supuesto se satisface” (Wooldridge, 2010, p.
48).
Para verificar en Python el cumplimiento del supuesto de variabilidad en las variables, se puede emplear,
por ejemplo, el método .apply(), utilizando como argumento la función std(...) de NumPy:
Nota: La función DataFrame(...) de pandas se usa en este código con fines puramente estéticos, para
obtener una presentación más agradable del resultado. La parte realmente importante, que genera el
contenido en el que se está interesado, es [Link]([Link], axis = 0).
Como se evidencia, todas las desviaciones estándar son diferentes de 0, por lo que hay variabilidad en
las variables (si la desviación estándar de alguna variable fuese 0, no se trataría de una variable, sino de una
constante).
In [7]: [Link]()
Si el usuario prefiere una representación un tanto más “rica” visualmente, puede recurrir a la función
heatmap(...) de Seaborn, la cual dibuja una matriz codificada por color. El argumento a pasar a la función,
para este caso específico, es la matriz de correlaciones, la cual, tal y como hemos apenas visto, es obtenida con
la aplicación del método .corr() sobre el DataFrame.
El código a emplear es el siguiente:
43
In [8]: fig, ax = [Link](figsize = (8, 6))
[Link]("Correlación simple\nentre variables explicativas",
fontsize = 18,
fontweight = "bold", x = 0.43)
[Link]([Link](), ax = ax)
plt.subplots_adjust(top = 0.85)
[Link](.9,-.02,
"Elaboración:",
fontsize = 13, fontweight = "bold",
ha = "right")
[Link](.9,-.07,
"Triana, F.\n(2019)",
fontsize = 12, ha = "right")
[Link]()
44
Nota: El lector debe notar que, dado que Seaborn es una librería basada en Matplotlib y que la función
heatmap(...) de Seaborn devuelve un Axes de una figura de Matplotlib, el código empleado en el blo-
que anterior es código común de esta última librería. Dado que Matplotlib desempeña un rol auxiliar en
este trabajo, no se dedicarán esfuerzos a examinar detalladamente la estructura de su código, por lo que
al lector interesado en ésta se sugiere consultar recursos adicionales.
Podemos observar que la escala de color usada por defecto por la función heatmap(...) de Seaborn no
parece ser la más adecuada. Es posible modificar esto incluyendo un valor específico para el parámetro cmap,
asignando la escala que nos parezca la más indicada. A modo de ejemplo, en el siguiente código se empleará la
escala ’RdYlGn’ (rojos, amarillos, verdes):
45
Si consideramos que nuestro gráfico aún no es lo suficientemente informativo, podemos agregar el coefi-
ciente de correlación simple correspondiente a cada una de las celdas. Para alcanzar tal objetivo, debemos
incluir el parámetro annot en la función heatmap(...) de Seaborn y asignarle el valor de True (annot es un
parámetro de tipo booleano). Así:
46
fontsize = 13, fontweight = "bold",
ha = "right")
[Link](.9,-.07,
"Triana, F.\n(2019)",
fontsize = 12, ha = "right")
[Link]()
Así, notamos que la correlación más fuerte se presenta entre las variables ’CM’ y ’FLR’, con un coeficiente
superior a 0.8 (en valor absoluto); sin embargo, tal magnitud del coeficiente no debe ser una preocupación
mayor, en cuanto el Supuesto #8 se refiere a colinealidad entre variables explicativas y la variable ’CM’ es la
variable dependiente, mientras que ’FLR’ es una variable explicativa.
El Supuesto #8 hace referencia a colinealidad exacta; esto es cuando una variable explicativa tiene una re-
lación lineal exacta con otra. Cuando se presenta colinealidad exacta, los coeficientes de regresión son indeter-
minados y los errores estándar son infinitos (Gujarati y Porter, 2010).
47
En nuestro caso específico, ninguna de las variables explicativas presenta una correlación muy fuerte con
alguna otra variable explicativa distinta a ella misma, por lo tanto, no se presenta multicolinealidad perfecta y
es posible estimar los parámetros.
Otra técnica útil para detectar multicolinealidad es por medio del Factor de Inflación de Varianza (VIF,
por sus siglas en inglés). “VIFs altos reflejan un incremento en las varianzas de los coeficientes de regresión
estimados debido a colinealidad entre variables predictoras, en comparación con las que se obtendrían cuando
las predictoras son ortogonales” (Murray, Nguyen, Lee, Remmenga, y Smith, 2012).
Así, en caso de que no haya multicolinealidad, o ésta sea leve, se espera que la magnitud de los VIFs sea
pequeña y entre más fuerte sea la colinealidad más alto será el VIF; Gujarati y Porter (2010), basados en el
trabajo de Kleinbaum, Kupper y Muller (1988), sugieren una regla práctica: “si el FIV [VIF] de una variable es
superior a 10 (esto sucede si R2j excede 0.90), se dice que esa variable es muy colineal” (p.340).
¿Cuáles son los VIFs de las variables del dataset con el que se está trabajando? Para obtener estos valores,
se hará uso de la función variance_inflation_factor(...) del módulo stats.outliers_influence
de la librería Statsmodels.
La función variance_inflation_factor(...) recibe dos argumentos: la matriz que contiene las va-
riables explicativas (exog) y un índice (exog_idx) que señala qué variable es a la que corresponde el VIF calculado.
Esta función calcula el VIF para la variable señalada; sin embargo, es de interés conocer el VIF de cada una de
las variables explicativas.
Es posible aplicar la función variance_inflation_factor(...) de forma “manual” tantas veces como
variables de interés haya; para nuestro caso, hay 3 variables explicativas, por lo que la función solo debe apli-
carse 3 veces, una operación que no es excesivamente tediosa. Sin embargo, si se tuviera un número apreciable
de variables explicativas, resultaría molesto repetir la operación en múltiples ocasiones, por lo que resulta útil
poder “automatizar” tal tarea. Una opción para lograr esto es utilizar un ciclo for:
Out[11]: VIF
FLR 1.711845
PGNP 1.078306
TFR 1.645150
Así, se puede observar que ninguna de las variables explicativas exhibe un VIF elevado, por lo que no se
sugiere que alguna de estas variables sea muy colineal. Este resultado está en línea con lo que se obtuvo en la
matriz de correlaciones, en la cual se encontró que ninguno de los coeficientes de correlación simple era muy
alto.
48
relativamente estable, con los puntos distribuidos simétricamente alrededor del cero, sin presentar patrones
específicos.
Dado que aún no se ha creado el modelo ni llevado a cabo la estimación, no se cuenta con una instancia de
resultados a emplear. Por lo tanto, lo primero que se hará es realizar la construcción y estimación correspon-
dientes:
In [12]: Y = data["CM"]
X = data[["PGNP", "FLR"]]
Modelo = [Link](Y, sm.add_constant(X))
Resultados = [Link]()
print([Link]())
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
[2] The condition number is large, 6.77e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Ya se cuenta con el modelo y sus resultados, pero aún no hemos asignado los residuos a un objeto, por lo que
llevaremos a cabo tal asignación, de modo que se simplifique el código empleado. Asignaremos los residuos,
los cuales son un atributo de la instancia de resultados (el usuario puede verificarlo con la función dir(...)),
al objeto Residuos, así:
49
Del mismo modo, es posible asignar los valores estimados, que se obtienen con la aplicación del método
.predict() a un objeto específico:
50
Como se observa, no parece que se presente un patrón específico y los puntos están simétricamente distri-
buidos alrededor del cero, por lo que el análisis visual sugiere que la forma funcional adoptada no es errónea.
Otro aspecto a considerar respecto a la estructura del modelo es que ésta se mantenga a lo largo de la esti-
mación; es decir, que no sea inestable. Para verificar la estabilidad en los parámetros, o no cambio estructural,
puede emplearse la prueba CUSUM basada en residuos de MCO, la cual toma como hipótesis nula no cambio
estructural.
La función empleada en Python para la prueba CUSUM basada en residuos de MCO es
breaks_olsresid(...) de [Link]; sus argumentos incluyen los residuos del modelo esti-
mado por MCO y el parámetro ddof, con valor por defecto igual a 0 y referente al número de grados de libertad
empleados en la corrección de la varianza del error.
Ahora, emplearemos el objeto Residuos como argumento de la función breaks_olsresid(...); los
resultados obtenidos con la ejecución de la función serán almacenados en un objeto al que denominaremos
ResultadosTest, así:
51
(0.5191974341185455, 0.9503227705917948, [(1, 1.63), (5, 1.36), (10, 1.22)])
Podemos ver que el resultado que obtenemos de nuestro código es tan solo una tupla Python (...) que
contiene un conjunto de datos; esto se debe a que [Link] no tiene integrado un modo de presentación
más “elaborado” del resultado.
Para mejorar la presentación de los resultados del Test se puede hacer uso de una Series de pandas, de
modo que podamos organizar adecuadamente la información disponible. Procederemos a crear dicha Series
asignándole como valores los componentes del resultado del test y como index una lista Python [...] con sus
respectivos nombres. Así:
Asumiendo un nivel de significancia de 5 %, como p − value > α entonces no se rechaza la hipótesis nula
y se concluye que no se presenta cambio estructural (a un nivel de confianza de 95 %).
Ya se han evaluado, muy brevemente, algunos de los supuestos del modelo clásico de regresión lineal sobre
la estructura del modelo. Ahora, se examinará el cumplimiento de los supuestos referentes al término de error.
In [18]: [Link](Residuos)
Out[18]: 1.254552017826427e-13
El argumento que recibe la función mean(...) de NumPy es el vector de residuos, el cual, conveniente-
mente, hemos asignado de forma previa al objeto Residuos. El resultado obtenido es prácticamente cero, pues
52
se trata de un número, expresado en notación científica, que es muy pequeño; esto sugiere que el Supuesto #3
se cumple para este ejemplo.
Si bien hemos utilizado una función, puede resultar útil señalar al lector que también existe un método,
con el que es posible calcular el promedio, el cual, como es de esperarse, es .mean():
In [19]: [Link]()
Out[19]: 1.254552017826427e-13
53
Como se puede observar, los residuos no exhiben un patrón particular y la amplitud de sus cambios no pare-
ce variar significativamente, por lo que, a primera vista, no hay evidencia de heteroscedasticidad. Sin embargo,
tal conclusión puede resultar apresurada.
Un análisis similar puede llevarse a cabo empleando los valores estimados y los residuos al cuadrado, lo
que permite controlar por la naturaleza de los residuos (positivos o negativos) al basarse en su magnitud y no
en su signo. Se espera que no se presente ningún patrón particular, de modo que la varianza de los residuos sea
constante y estos no dependan de la magnitud de la variable regresada.
Para ejecutar este análisis utilizaremos el método .scatter(...) empleando como argumentos el vec-
tor de valores estimados y el de residuos al cuadrado. Ya contamos con el primero, pero solo se tiene infor-
mación de los residuos y no de los residuos cuadrados; sin embargo, basta con aplicar el operador ** para
obtener dichas magnitudes, ya que éste es el que permite en Python realizar potenciación siguiendo la estruc-
tura base**exponente:
54
fontsize = 13, fontweight = "bold",
ha = "right")
[Link](.9,-.08,
"Triana, F.\n(2019)",
fontsize = 12, ha = "right")
plt.subplots_adjust(top = 0.85)
[Link]()
De acuerdo al resultado obtenido, no parece que los residuos sean heteroscedásticos, aunque hay una li-
gera impresión de que van creciendo con la magnitud de la variable dependiente; para llegar a una conclusión
“válida”, el análisis visual no es suficiente: se requiere de pruebas estadísticas formales.
Una de las pruebas estadísticas generalmente empleadas para verificar la presencia de heteroscedasticidad
es el Test Breusch-Pagan de Multiplicadores de Lagrange (Greene, 2003), el cual, básicamente, toma como
hipótesis nula homoscedasticidad.
La función empleada en Python para el Test Breusch-Pagan de Multiplicadores de Lagrange para heteros-
cedasticidad es het_breuschpagan(...) de [Link]; sus argumentos incluyen los residuos del mo-
delo y el conjunto de variables que pueden causar la heteroscedasticidad.
Ya se cuenta con un objeto que contiene los residuos, sin embargo, aún no se han asignado las variables
explicativas del modelo a un objeto particular. Como las variables explicativas consisten en un atributo de la
instancia de resultados del modelo, podemos asignarlas a un objeto aplicando dicho atributo, así:
Ya que contamos con los objetos a emplear como argumentos, podemos hacer uso de la función
het_breuschpagan(...) de [Link]. Siguiendo la estructura que hemos empleado para las pruebas
de la sección anterior, el código a ejecutar es el siguiente:
55
In [23]: ResultadosTest = sms.het_breuschpagan(Residuos, Explicativas)
Nombres = ["Estadístico LM", "p-value del estadístico LM", "Estadístico F",
"p-value del estadístico F"]
[Link](ResultadosTest, index = Nombres)
Asumiendo un nivel de significancia de 5 %, como p − value > α entonces no se rechaza la hipótesis nula
y se concluye que se presenta homoscedasticidad en el término de error (trabajando con un α de 0.05).
Nota 1: Este test equivale al generado por la función bptest(...) de R (paquete lmtest).
Nota 2: Este test equivale al generado por el código estat hettest, rhs fstat de Stata.
Otra de las pruebas estadísticas comúnmente empleadas para verificar la presencia de heteroscedasticidad
es el Test White de Multiplicadores de Lagrange, el cual hace uso de regresiones auxiliares, toma, esencialmente,
como hipótesis nula homoscedasticidad, y es capaz de identificar formas más generales de heteroscedasticidad
que el Test de Breusch-Pagan (Verbeek, 2004).
La función empleada en Python para el Test White de Multiplicadores de Lagrange para heteroscedasticidad
es het_white(...) de [Link]; sus argumentos incluyen los residuos del modelo y el conjunto de
variables a emplear en las regresiones auxiliares.
Los argumentos empleados en la función het_white(...) de [Link] corresponden a los mismos
utilizados en la función het_breuschpagan(...) de dicho módulo y, dado que estos han sido asignados
previamente a objetos específicos, basta con reutilizar tales objetos.
Siguiendo la estructura que se ha empleado para las demás pruebas estadísticas realizadas, el código a
ejecutar es el siguiente:
Los resultados del Test White de Multiplicadores de Lagrange para heteroscedasticidad reafirman la con-
clusión del Test Breusch-Pagan realizado con anterioridad. Asumiendo un nivel de significancia de 5 %, como
p − value > α entonces no se rechaza la hipótesis nula y se concluye que se presenta homoscedasticidad en
el término de error (a un nivel de confianza de 95 %).
Nota: Este test equivale al generado por el código estat imtest, white de Stata.
56
No autocorrelación, o correlación serial, entre las perturbaciones (Supuesto #5)
El supuesto de no autocorrelación entre las perturbaciones indica que estas no siguen patrones sistemáticos
y no están correlacionadas entre sí; tal supuesto puede justificarse para el caso de datos transversales, pero
tiende a incumplirse cuando se trabaja con series de tiempo, en cuanto las observaciones sucesivas usualmente
están fuertemente correlacionadas (Véase la Sección 3.2, del Capítulo 3, Modelos de regresión con dos variables:
problema de estimación, de Gujarati y Porter (2010)).
Este supuesto normalmente se verifica con pruebas estadísticas, algunas de las cuales se tratan a continua-
ción. Una prueba estadística generalmente empleada para verificar la presencia de autocorrelación de orden 1
es el Test Durbin-Watson, el cual toma, esencialmente, como hipótesis nula no autocorrelación (de orden 1) y
genera un estadístico cuyo valor se contrasta con los valores críticos correspondientes (Chatterjee y Simonoff,
2013) para establecer una conclusión.
La función empleada en Python para el Test Durbin-Watson es durbin_watson(...) de [Link];
sus argumentos incluyen los residuos del modelo.
Siguiendo la misma estructura de las otras pruebas estadísticas realizadas, y teniendo en cuenta que los
residuos ya se asignaron a un objeto específico, el código a ejecutar es el siguiente:
El valor obtenido es 2.1861, entonces, ¿se presenta autocorrelación de orden 1? El lector debe notar que, a
diferencia de otras pruebas estadísticas, en el caso del Test Durbin-Watson no se reporta un p-value con el cual
decidir si se rechaza o no la hipótesis nula, entonces, ¿cómo se determina la conclusión?
El estadístico Durbin-Watson es contrastado con valores críticos a partir de los cuales se establecen zonas
de autocorrelación positiva, no autocorrelación, indeterminación y autocorrelación negativa; dependiendo de
la zona en la que se ubique el valor del estadístico, se obtiene una conclusión específica.
Para este caso en concreto, asumiendo un nivel de significancia de 5 %, los valores de d L y dU para 2 variables
explicativas (excluyendo el término constante) y 65 observaciones son, respectivamente, 1.536 y 1.662. El
valor de 2.1861 pertenece al intervalo (dU , 4 − dU ), el cual corresponde a la zona de No Autocorrelación, por
lo que se concluye que no se presenta autocorrelación de primer orden (trabajando con un α de 0.05).
Nota: Este test equivale al generado por la función dwtest(...) de R (paquete lmtest). El usuario debe
tener en cuenta que, en la función de R, a diferencia de la de Python, el argumento a emplear no son los
residuos sino el propio modelo.
El Test Durbin-Watson es una prueba estadística que permite identificar autocorrelación de primer orden,
sin embargo, si se quiere examinar correlación serial de orden superior ésta ya no puede emplearse y gene-
ralmente se recurre al Test Breusch-Godfrey, el cual toma, esencialmente, no autocorrelación como hipótesis
nula y se basa en la ejecución de una regresión auxiliar en la que se incluyen rezagos del término de error como
variables explicativas (Kleiber y Zeileis, 2008).
La función empleada en Python para el Test Breusch-Godfrey de Multiplicadores de Lagrange para autoco-
rrelación de los residuos es acorr_breusch_godfrey(...) de [Link]; sus argumentos incluyen la
instancia de resultados (¡no los residuos!) del modelo y el número de rezagos, nlags, a incluir en la regresión
auxiliar.
Siguiendo la misma estructura de las otras pruebas estadísticas realizadas, y teniendo en cuenta que la
instancia de resultados del modelo ya se asignó a un objeto específico, el código a ejecutar es el siguiente:
57
In [26]: ResultadosTest = sms.acorr_breusch_godfrey(Resultados, nlags = 1)
Nombres = ["Estadístico LM", "p-value del estadístico LM", "Estadístico F",
"p-value del estadístico F"]
[Link](ResultadosTest, index = Nombres)
Los resultados del Test Breusch-Godfrey de Multiplicadores de Lagrange para autocorrelación de los resi-
duos reafirman la conclusión del Test Durbin-Watson realizado con anterioridad. Asumiendo un nivel de signi-
ficancia de 5 %, como p − value > α entonces no se rechaza la hipótesis nula y se concluye que no se presenta
correlación serial de primer orden en el término de error (trabajando con un nivel de significancia de 5 %).
Ya se ha indicado que el Test Breusch-Godfrey es una prueba general de autocorrelación, en cuanto permite
identificar correlación serial de orden superior a 1. En el código anterior hemos asignado el valor de 1 al pará-
metro nlags de la función acorr_breusch_godfrey(...) de [Link], por lo que hemos examinado
la presencia de autocorrelación de primer orden, al igual que en el Test Durbin-Watson. Ahora, procederemos
a examinar la presencia de autocorrelación hasta de orden 4, por lo que nlags tomará el valor de 4. Basta con
reutilizar la última celda de código que hemos empleado y modificar el valor del parámetro nlags, así:
Nota: Este test equivale al generado por la función bgtest(...) de R (paquete lmtest).
58
In [28]: ResultadosTest = sms.acorr_ljungbox(Residuos, lags = 1, boxpierce = True)
Nombres = ["Estadístico LB", "p-value del estadístico LB", "Estadístico BP",
"p-value del estadístico BP"]
[Link](ResultadosTest, index = Nombres)
Tanto el p-value correspondiente al estadístico del Test Ljung-Box como el correspondiente al Test Box-
Pierce son superiores a 0.05, por lo que no se rechaza la hipótesis nula y se concluye que no se presenta corre-
lación serial de primer orden (observe que lags = 1).
Al igual que en el Test Breusch-Godfrey, el Test Ljung-Box y el Test Box-Pierce pueden emplearse para
examinar la existencia de correlación serial hasta de orden p. Para examinar hasta p orden de autocorrelación,
basta con asignar el valor de p al parámetro lags de la función acorr_ljungbox(...) de [Link].
Siguiendo la misma estructura de las otras pruebas estadísticas realizadas, y teniendo en cuenta que los
residuos ya se asignaron a un objeto específico, el código a ejecutar para evaluar correlación serial hasta de
orden 2 (por ejemplo) es el siguiente:
El lector debe notar que la función acorr_ljungbox(...) de [Link] genera p valores para el
estadístico y p-value correspondientes, los cuales presenta dentro de una lista Python “[...]”, en cuanto ejecuta
la prueba respectiva para cada uno de los rezagos incluidos.
Nota: Este test equivale al generado por la función [Link](...) de R (paquete stats). A diferencia
de Python, en R la prueba que se ejecuta por defecto es la de Box-Pierce, por lo que si se quiere llevar a
cabo el Test Ljung-Box, tal preferencia debe señalarse explícitamente en el parámetro type de la función.
Supuesto de Normalidad
Uno de los supuestos sobre el término de error es que tiene distribución de probabilidad normal. Es posible
que el lector cuidadoso haya notado que este supuesto no se incluyó en el listado de supuestos del modelo
clásico de regresión lineal, presentado en la primera parte de este capítulo. Esto se debe a que el supuesto de
normalidad en la distribución del término de error no hace parte del modelo clásico de regresión lineal, por
lo que no es necesario para garantizar las propiedades de los estimadores de MCO de linealidad, insesgadez y
varianza mínima. Entonces, ¿por qué se menciona?
59
El supuesto de normalidad en la distribución del término de error hace parte del modelo clásico de regresión
lineal normal, y ve justificada su existencia por la necesidad de llevar a cabo inferencia estadística. Mientras los
demás supuestos presentados a inicio de este capítulo son requeridos para que los estimadores de MCO sean
MELI (o BLUE en inglés), tales supuestos no son suficientes para garantizar que la inferencia estadística sea
válida.
Debido a que el método de MCO no hace ninguna suposición respecto de la naturaleza probabilís-
tica de ui , resulta de poca ayuda para el propósito de hacer inferencias sobre la FRP mediante la
FRM, a pesar del teorema de Gauss-Markov. Este vacío puede llenarse si se supone que las u siguen
una determinada distribución de probabilidad. Por razones que mencionaremos en seguida, en el
contexto de regresión se supone, por lo general, que las u tienen la distribución de probabilidad
normal. (Gujarati y Porter, 2010, p. 98)
El supuesto indica que el término de error tiene distribución normal con media cero y varianza constante
y que, además, la covarianza entre ui y u j (i 6= j) es cero. Así, el supuesto puede resumirse en que ui ∼
N ID (0, σ2 ), en otras palabras: el término de error es normal e independientemente distribuido.
Una de las pruebas estadísticas más populares empleadas para verificar el cumplimiento de normalidad en
la distribución del término de error es el Test Jarque-Bera, el cual se basa en la asimetría y la curtosis (Würtz y
Katzgraber, 2009) y toma como hipótesis nula, básicamente, distribución normal.
La función empleada en Python para el Test Jarque-Bera para normalidad es jarque_bera(...) de stats-
[Link]; sus argumentos incluyen los residuos del modelo.
Siguiendo la misma estructura de las otras pruebas estadísticas realizadas, y teniendo en cuenta que los
residuos ya se asignaron a un objeto específico, el código a ejecutar para evaluar normalidad en la distribución
del término de error es el siguiente:
Asumiendo un nivel de significancia de 5 %, como p − value > α entonces no se rechaza la hipótesis nula
y se concluye que se presenta normalidad en la distribución del término de error (trabajando con un α de 0.05).
Nota: Este test equivale al generado por la función [Link](...) de R (paquete tseries). El
usuario debe tener en cuenta que el argumento empleado en la función de R, al igual que la de Python,
corresponde a los residuos del modelo.
Otras pruebas ampliamente utilizadas para evaluar normalidad son la de Anderson-Darling, la cual perte-
nece a la clase cuadrática de los estadísticos basados en Función de Distribución Empírica, y la de Kolmogorov-
Smirnov, también construida a partir de dicha función (Razali y Wah, 2011). La primera prueba es implementada
en Python con la función normal_ad(...) de [Link]:
60
C:\Users\FCE\Anaconda3\lib\site-packages\statsmodels\stats\_adnorm.py:66:
FutureWarning: Using a non-tuple sequence for multidimensional indexing
is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`.
In the future this will be interpreted as an array index, `arr[[Link](seq)]`,
which will result either in an error or a different result.
S = [Link]((2*i[sl1]-1.0)/N*([Link](z)+[Link](1-z[sl2])), axis=axis)
Reforzando la conclusión del Test Jarque Bera, el resultado de Anderson-Darling sugiere que la distribución
de los residuos es normal.
Respecto a la otra prueba, la de Kolmogorov-Smirnov, esta es implementada con la función
kstest_normal(...) del módulo diagnostic de [Link]. La estructura a emplear para presen-
tar los resultados es la misma utilizada en las demás pruebas estadísticas que se han realizado:
Al igual que con Jarque-Bera y Anderson-Darling, los resultados del test de Kolmogorov-Smirnov indican
que no se rechaza la hipótesis nula y, por tanto, no se rechaza que la distribución del término de error sea
normal.
Las dos últimas pruebas realizadas se basan en la Función de Distribución Empírica, la cual, de acuerdo a
Vaart (1998) es el estimador natural de la distribución subyacente, F, si esta no es conocida. Para visualizar
la distribución empírica en Python podemos hacer uso de funciones de Matplotlib y para obtener los valores
correspondientes a dicha función basta con hacer unas cuantas operaciones sencillas.
Vaart (1998) define, para una muestra aleatoria X1 , ..., Xn , la Función de Distribución Empírica de la si-
guiente manera:
1 n
n i∑
Fn (t) = 1 { Xi ≤ t }
=1
A partir de esta definición es muy simple determinar los puntos con los que se calculará la Función de
Distribución Empírica. En primer lugar, debemos organizar los valores de la variable de interés de menor a
mayor, para lo que se utiliza la función sort(...) de NumPy:
In [33]: x = [Link](Residuos)
A cada uno de estos valores x corresponde una “probabilidad” de ser menor o igual; es decir, si tenemos el
conjunto {1,2,2,3,5,7,8,9,10,10} la probabilidad de que un elemento sea menor o igual a 2 es 30 % (pues tres
de los diez elementos son menores o iguales a dos), la probabilidad de que un elemento sea menor o igual a 7
es 60 % (pues seis de los diez elementos son menores o iguales a siete). Esta lógica es muy fácil de implementar
en el código: tan solo se debe determinar el número de elementos, n, y dividir una secuencia 1, ..., n entre este
número, de modo que se obtengan las probabilidades correspondientes:
61
In [34]: n = [Link]
y = [Link](1, n+1) / n
Ya se tienen los pares ( x, y) para graficar la Función de Distribución Empírica, pero ésta solo nos resulta
útil si podemos compararla con algún referente: una distribución teórica. La distribución normal se construye
con dos parámetros: la media y la desviación estándar. Si se generan datos a partir de una distribución normal
con media igual al promedio de los residuos y desviación estándar igual a la de estos, dichos datos sirven para
graficar la distribución téorica (pues se sabe que efectivamente se trata de una distribución normal), la cual
puede usarse para comparar con la distribución empírica; si el comportamiento es muy similar, visualmente se
sugiere que la distribución de los residuos es normal.
Veamos lo anteriormente descrito en la práctica. Para generar valores aleatorios provenientes de una dis-
tribución normal se puede usar la función normal(...) del módulo random de NumPy. A esta función se debe
especificar la media (loc), la desviación estándar (scale) y el número de valores a generar (size):
A partir de estos residuos teóricos puede definirse la función de distribución “teórica”, para lo que basta
emplear las mismas operaciones usadas en la distribución empírica:
62
plt.subplots_adjust(top = 0.85)
[Link]()
El lector puede observar que la distribución de los residuos no dista mucho de la apariencia que corresponde
a la distribución normal. Los puntos correspondientes a la función de distribución empírica (los azules) se sitúan
muy cerca de la línea de referencia (la distribución teórica).
Finalmente, otro método gráfico ampliamente utilizado, seguramente más común que el de la distribución
empírica, para evaluar la normalidad, es el gráfico de cuantil-cuantil, o Q-Q Plot. Este gráfico puede construirse
con la función qqplot(...) de Statsmodels, pasando como argumento el objeto que contiene los datos a partir
de los cuales se quieren establecer los cuantiles. Por defecto esta función trabaja con la distribución normal,
razón por la cual no debe introducirse alguna modificación:
63
color = "darkblue")
ax.set_xlabel("Cuantiles teóricos", fontsize = 14)
ax.set_ylabel("Cuantiles muestrales", fontsize = 14)
[Link](.9,-.02,
"Elaboración:",
fontsize = 13, fontweight = "bold",
ha = "right")
[Link](.9,-.08,
"Triana, F.\n(2019)",
fontsize = 12, ha = "right")
[Link]()
Como se observa, la mayoría de puntos se sitúa muy cerca de la línea roja de referencia, indicador de que
la distribución de los datos no dista mucho de la teórica (la normal). Este resultado concuerda con el obtenido
en el análisis de la distribución empírica y con los resultados de las tres pruebas estadísticas empleadas para
verificar el cumplimiento del supuesto de normalidad.
64
En este capítulo se abordó muy brevemente, y sin examinar detalladamente los fundamentos teóricos aso-
ciados, el cumplimiento de algunos de los supuestos del modelo clásico de regresión lineal (normal, si consi-
deramos el último supuesto planteado) usando herramientas de librerías Python especializadas. El modelo al
que corresponden tales supuestos se ilustró, por medio de un ejemplo particular, en el Capítulo 3, en el cual se
dio una ligera aproximación (práctica) al tema fundamental que generalmente se aborda en los cursos de nivel
introductorio de Econometría. Ahora, se aumentará ligeramente la dificultad procediendo a examinar un tema
un tanto más complejo (en un sentido más teórico que práctico) al que se dedicará el próximo capítulo.
65
Capítulo 5
Uno de los supuestos (el #2 específicamente) del modelo clásico de regresión lineal planteados a inicios del
capítulo previo, establece que la covarianza entre las variables explicativas y el término de error es cero, es decir,
que toda x (siguiendo la notación empleada con anterioridad) es independiente del término de error. Cuando
las variables explicativas cumplen esta condición, son denominadas exógenas y cuando violan tal supuesto se
consideran endógenas.
Cuando alguna(s) de las variables explicativas del modelo es (son) endógena(s) los estimadores obtenidos
por el método de Mínimos Cuadrados Ordinarios pierden su calidad de MELI. Para solucionar tal inconveniente
se recurre al método de Variables Instrumentales, en el cual se emplean variables no incluidas en la ecuación
original para estimar las variables explicativas endógenas incluidas en ésta y tratar su endogeneidad.
Cuando alguna(s) variable(s) explicativa(s) y el término de error se correlacionan, se violan condiciones
de los supuestos del modelo clásico de regresión lineal, por lo que los estimadores de MCO pierden algunas
de sus características atractivas. Para solucionar el inconveniente generado por la endogeneidad de algún(os)
covariante(s) se requiere de información adicional, la cual se obtiene de una variable observable exógena no
incluida.
Tomando un modelo de regresión lineal de la forma
y1i = β 0 + β 1 xi + β 2 y2i + ui
en donde y1 es la variable dependiente, x es una variable explicativa exógena (no se correlaciona con el
término de error) y y2 es una variable explicativa endógena (se correlaciona con el término de error), si se tiene
una variable z que cumple con las condiciones siguientes:
1. cov(z, u) = 0
2. cov(z, y2 ) 6= 0
66
variable(s) explicativa(s); sin embargo, aún no se ha señalado el proceso específico por medio del cual se hace
uso de los instrumentos.
Uno de los métodos empleados para obtener estimadores de Variables Instrumentales es el de Mínimos
Cuadrados en 2 Etapas (MC2E), el cual recibe tal denominación dado que el proceso que comprende está
estructurado en dos “etapas” bien definidas en las que se llevan a cabo regresiones específicas con objetivos
concretos.
Tomando el modelo de regresión lineal señalado anteriormente, y siguiendo a Wooldridge (2010), la pri-
mera etapa del método MC2E consiste en realizar la regresión ŷ2i = π̂0 + π̂1 xi + π̂2 zi , de donde se ob-
tienen los valores ajustados ŷ2 . La segunda etapa es una regresión por Mínimos Cuadrados Ordinarios de
y1i = β 0 + β 1 xi + β 2 ŷ2i + ui (nótese que en esta regresión se emplea los valores de ŷ2 en vez de los va-
lores de y2 ). Así, de lo que se trata, básicamente, es de estimar la variable explicativa endógena a partir de
variables exógenas y, posteriormente, emplear la variable estimada como variable explicativa, junto a las de-
más exógenas del modelo original (no se incluyen los instrumentos), para estimar la variable dependiente.
Ya se ha dado una muy mínima descripción de Variables Instrumentales y el método de Mínimos Cuadrados
en 2 Etapas, por lo que el lector debe tener una idea general del tema que ahora se pretende abordar de manera
práctica; labor a la cual, a continuación, se da inicio efectivo con Python.
El dataset que se utilizará para ilustrar el uso de variables instrumentales y el método MC2E es el empleado
en el Ejemplo 15.5, Rendimientos de la educación para la mujer trabajadora, de Wooldridge (2010). El dataset se
ha obtenido desde Stata, con el comando bcuse, y ha sido exportado como archivo .csv para su uso en Python.
Se invita al lector a consultar el ejemplo, de modo que pueda verificar los resultados obtenidos.
El lector ya debe saber que el anterior bloque de código permite cargar las herramientas necesarias para
llevar a cabo las acciones requeridas en el análisis econométrico planteado; en caso de que no tenga claridad
al respecto, se le pide consultar la sección ’Preparación del entorno’ del segundo capítulo, en donde se brinda
una explicación sobre tal cuestión.
El código de esta celda es exactamente igual al empleado en la preparación del entorno del Capítulo 2,
exceptuando por lo siguiente: 1) no se importa la librería Seaborn y, 2) más importante aún, se importa la función
IV2SLS(...) del módulo iv de la librería linearmodels; tal función es fundamental para el tema a abordar en
este capítulo y su importación es obligatoria.
El usuario debe ser muy cuidadoso al ejecutar la celda de código anterior, pues si lo hace inmediatamente,
sin tener en cuenta la información que se presenta a continuación, el proceso no será completamente exitoso;
esto, por la razón que enseguida se presenta. La librería linearmodels no viene integrada por defecto en la Distri-
bución Anaconda (la cual es la que estamos empleando), por lo cual es necesario instalarla por nuestra propia
cuenta para poder hacer uso de las herramientas que brinda.
Para instalar una librería basta con dirigirse al Símbolo del Sistema (’cmd’); estando ahí, tan solo se debe
escribir conda install y, a continuación, el nombre de la librería en la que se está interesado. En nuestro
67
caso, el código a emplear es conda install linearmodels. De forma alternativa puede emplearse !conda
install linearmodels directamente en el notebook.
Nota: El código anterior hace uso de conda, el gestor de paquetes propio de la distribución Anaconda, el
cual recurre al repositorio Anaconda; sin embargo, es posible que el proceso de instalación de la librería
a través de conda no resulte exitoso, en cuyo caso se puede utilizar el sistema de gestión de paquetes
estándar de Python, pip, el cual trabaja directamente con Python Package Index (PyPi). El código para la
instalación de linearmodels a través de pip es pip install linearmodels en ’cmd’ o !pip install
linearmodels en el notebook; ejecutando dicho código debería llevarse a cabo de forma exitosa la ins-
talación de la librería linearmodels.
La importación de linearmodels es fundamental para nuestra labor, pues es donde están contenidas las fun-
ciones que emplearemos para hacer uso de variables instrumentales y aplicar el método de Mínimos Cuadra-
dos en 2 Etapas; sin esta librería resultaría excesivamente complejo el proceso de estimación y la realización
de las labores requeridas para llevar a cabo un análisis econométrico “aceptable”. El usuario debe verificar la
instalación previa de linearmodels (en caso de no estar instalada, llevar a cabo el proceso de instalación corres-
pondiente, sea con conda o con pip) y, en caso de reutilizar el bloque de código de preparación del entorno
del segundo capítulo, debe ser cuidadoso y asegurarse de incluir correctamente la línea de importación de la
función IV2SLS(...) del módulo iv de linearmodels.
Una vez ejecutada la celda de código de la preparación del entorno, sin la generación de algún tipo de error,
puede continuarse con entera tranquilidad a la fase siguiente, la cual tiene una importancia fundamental en
cuanto es la que permite contar con la “materia prima” para llevar a cabo el análisis: los datos.
Al lector que no entienda esta línea de código, se le sugiere dirigirse al Capítulo 2, en donde encontrará una
explicación que seguramente aclarará sus dudas.
Una vez ejecutada la línea de código de esta celda, el dataset debería haber sido importado correctamente.
Para comprobarlo, recurrimos al método .head(), aplicado al objeto que contiene el dataset (en este caso le
hemos asignado el nombre data):
In [3]: [Link]()
Out[3]: inlf hours kidslt6 kidsge6 age educ wage repwage hushrs husage \
0 1 1610 1 0 32 12 3.3540 2.65 2708 34
1 1 1656 0 2 30 12 1.3889 2.65 2310 30
2 1 1980 1 3 35 12 4.5455 4.04 3072 40
3 1 456 0 3 34 12 1.0965 3.25 1920 53
4 1 1568 1 2 31 14 4.5918 3.60 2000 32
68
0 … 16310 0.7215 12 7 5.0 0 14 10.910060
1 … 21800 0.6615 7 7 11.0 1 5 19.499981
2 … 21040 0.6915 12 7 5.0 0 15 12.039910
3 … 7300 0.7815 7 7 5.0 0 6 6.799996
4 … 27300 0.6215 12 14 9.5 1 7 20.100060
lwage expersq
0 1.210154 196
1 0.328512 25
2 1.514138 225
3 0.092123 36
4 1.524272 49
[5 rows x 22 columns]
Podemos observar que el proceso de importación ha sido exitoso y la información del dataset ha sido al-
macenada correctamente en un DataFrame de pandas. Ahora, podemos continuar tranquilamente con el desa-
rrollo de nuestro análisis, en cuanto hemos concluido satisfactoriamente el primer paso. Procederemos a usar
variables instrumentales empleando el método de Mínimos Cuadrados en 2 Etapas, sin embargo, el usuario
debe tener en cuenta la información que se le brinda a continuación.
El dataset que acabamos de importar contiene algunas observaciones faltantes (esto se debe a la fuente
original de la información y no a una preferencia del autor de esta obra, aunque veremos que tal situación nos
resultará algo enriquecedora al mejorar nuestra capacidad exploratoria y de manejo de datos). Estos valores
faltantes pueden generar problemas posteriormente, por lo que procederemos a identificarlos y darles un tra-
tamiento apropiado.
Para examinar cuántas filas presentan algún campo con un valor faltante y conocer su ubicación, emplea-
remos los métodos .isna() y .sum(); el primero nos permite identificar valores faltantes (por medio de un
criterio booleano) y el segundo permite ejecutar una suma que, teniendo en cuenta el tipo de datos generados
por el primer método, permite conocer la cantidad de los mismos.
In [4]: [Link]().sum()
Out[4]: inlf 0
hours 0
kidslt6 0
kidsge6 0
age 0
educ 0
wage 325
repwage 0
hushrs 0
husage 0
huseduc 0
huswage 0
faminc 0
mtr 0
motheduc 0
fatheduc 0
69
unem 0
city 0
exper 0
nwifeinc 0
lwage 325
expersq 0
dtype: int64
Se puede observar que las únicas variables en las que se presentan valores faltantes son ’wage’ y ’lwage’,
con 325 datos clasificados como Na para cada una. Esta situación puede dar lugar a ciertos problemas, por lo
cual, tratando de minimizar la posibilidad de inconvenientes durante la fase de estimación, solo se empleará
los registros en los cuales existe un valor específico (no Na) para la variable ’lwage’ (la cual tiene el rol de
variable dependiente en el modelo). ¿Cómo podremos conseguir tal objetivo, seleccionar únicamente registros
con valores no faltantes?
Para obtener un DataFrame que solamente contenga registros en los que no haya valores faltantes para
la variable ’lwage’ se hará uso de la selección por medio de corchetes “[ ]” y del método .notna() (aplicado
sobre la variable de referencia); esto, con el propósito de tener en cuenta únicamente registros que cumplan la
condición señalada. El resultado obtenido se asignará al objeto data (que ya existe), con el objetivo de contar
con un único dataset y evitar posibles confusiones. El código a ejecutar es el siguiente:
In [5]: data = data[data["lwage"].notna()]
Ahora, procedemos a verificar que hemos conseguido lo que nos hemos propuesto y en nuestro DataFrame
no se presentan valores faltantes; de nuevo, emplearemos los métodos .isna() y .sum():
In [6]: [Link]().sum()
Out[6]: inlf 0
hours 0
kidslt6 0
kidsge6 0
age 0
educ 0
wage 0
repwage 0
hushrs 0
husage 0
huseduc 0
huswage 0
faminc 0
mtr 0
motheduc 0
fatheduc 0
unem 0
city 0
exper 0
nwifeinc 0
lwage 0
expersq 0
dtype: int64
70
Observamos que todas las variables contenidas en el DataFrame carecen de valores faltantes, por lo que
ningún registro está incompleto. Este DataFrame, que tiene el nombre de data, es el que emplearemos para
trabajar con variables instrumentales y aplicar el método de Mínimos Cuadrados en 2 Etapas, el cual, aborda-
remos a continuación desde dos enfoques: 1) proceso manual y 2) proceso automático.
Proceso manual
El método de Mínimos Cuadrados en 2 Etapas, tal y como su nombre lo indica, está formado por dos etapas
de estimación claramente definidas: la primera consiste en regresar la(s) variable(s) explicativa(s) endógena(s)
contra el(los) instrumentos y la(s) variable(s) explicativa(s) exógena(s); la segunda etapa consiste en emplear
la variable estimada en la primera etapa en lugar de la variable explicativa endógena en el modelo original. Al
lector que tal descripción le parezca confusa, se sugiere remitirse a la primera parte de este capítulo o consultar
el Capítulo 15, Estimación con variables instrumentales y mínimos cuadrados en dos etapas, de Wooldridge (2010),
donde este tema es tratado en detalle.
El proceso “manual” recibe tal denominación en cuanto debemos llevar a cabo las regresiones de la primera
y la segunda etapa por nuestra propia cuenta, especificando el conjunto de variables a emplear en cada una
de estas. Tal procedimiento no es recomendable y la mayoría de programas y lenguajes que soportan análisis
estadístico (Python no es la excepción) cuenta con funciones que realizan la ejecución “automáticamente”,
evitando que el usuario deba llevar a cabo las tareas por su cuenta.
A pesar de lo señalado en el párrafo anterior, dado el enfoque práctico de esta guía y la utilidad que puede
tener para el usuario el llevar a cabo el proceso paso por paso en la comprensión del tema tratado, procederemos
a aplicar el método MC2E ’manualmente’ (el proceso automático se aborda en la siguiente sección, por lo que
al lector que no esté interesado en ejecutar ambas etapas por su propia cuenta se le sugiere consultar dicha
sección directamente, sin dedicar su tiempo al estudio de la presente).
Lo primero que debemos hacer para llevar a cabo correctamente las regresiones de la primera etapa y la
segunda etapa es agrupar las variables disponibles de forma apropiada, en cuanto cada una de estas desem-
peña un rol específico. Así, procederemos a asignar la variable dependiente a un objeto específico, la variable
explicativa endógena a otro objeto en particular y, del mismo modo, las variables explicativas exógenas y los
instrumentos a otros objetos concretos. El código a ejecutar es el siguiente:
Primera etapa
Durante la primera etapa se regresa la variable explicativa endógena contra las variables explicativas exóge-
nas y los instrumentos. El lector puede observar que las variables explicativas exógenas y los instrumentos han
sido asignados a objetos distintos, por lo que, con el propósito de simplificar el código, ahora se asignarán
conjuntamente a un único objeto. Para conseguir tal objetivo, basta con emplear la función concat(...) de
pandas, usando como argumento una lista Python “[...]” en la que se incluyan los DataFrames que contienen
las variables explicativas exógenas y los instrumentos, y asignando el valor de 1 al parámetro axis (para conca-
tenar los DataFrames “horizontalmente”). El resultado generado por la aplicación de la función sobre la lista
recibirá el nombre (algo extenso, pero muy informativo) de RegresorasPrimeraEtapa, así:
71
Verificamos que el proceso se ha llevado a cabo exitosamente, recurriendo al método .head():
In [9]: [Link]()
Ahora, dado que nuestras variables están almacenadas de forma conveniente, construimos el modelo co-
rrespondiente a la primera etapa, tomando como variable dependiente la variable explicativa endógena (’educ’
en notación original, Y2 en nuestro entorno) y como variables explicativas las variables explicativas exóge-
nas (’exper’, ’expersq’ en notación original) y los instrumentos (’motheduc’, ’fatheduc’ en notación origi-
nal). Nótese que las variables explicativas exógenas y los instrumentos están contenidos en un único objeto
(RegresorasPrimeraEtapa).
El código a ejecutar el siguiente:
72
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
[2] The condition number is large, 1.55e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Al lector que ha seguido atentamente esta obra, la anterior celda de código no debería resultarle extraña,
pues no es más que una regresión por Mínimos Cuadrados Ordinarios; si el código presentado le resulta confuso,
se le sugiere consultar el tercer capítulo de este trabajo, en donde la regresión por MCO es el tema central.
Hemos empleado instrumentos como covariantes en una regresión, sin embargo, aún desconocemos si di-
chos instrumentos cumplen las condiciones para desempeñar satisfactoriamente tal papel; todavía no hemos
evaluado la “calidad” de las variables instrumentales utilizadas, por lo que procederemos a examinarla.
Para que las variables instrumentales empleadas cumplan la condición de relevancia (véase la primera parte
de este capítulo), los parámetros asociados a los instrumentos deben ser estadísticamente distintos de cero.
Como puede observarse en la instancia de resultados del modelo, tanto para la variable ’motheduc’ como para
la variable ’fatheduc’ se cumple dicha condición, pues los coeficientes asociados son estadísticamente signi-
ficativos, el p-value respectivo es inferior a 0.05 (nivel de significancia elegido) y el intervalo de confianza no
contiene el 0.
La conclusión del párrafo anterior es que los instrumentos cumplen la condición de relevancia, sin embargo,
para probar formalmente la validez de dicha conclusión puede emplearse una Prueba F para significancia de un
subconjunto de parámetros. La Prueba F se ejecuta en Python por medio del método .f_test(), aplicado sobre
la instancia de resultados del modelo y empleando como argumento las condiciones a evaluar (que pueden
expresarse como un ’string’).
En este caso, asignaremos la hipótesis a un objeto que se empleará como argumento en el método
.f_test(); esto, con el propósito de simplificar el código. Así:
El resultado generado por el método .f_test() es un objeto de tipo ContrastResults(), el cual cuenta
con sus propios atributos. Algunos de estos atributos corresponden al estadístico y al p-value, los cuales son
los valores de interés para nuestros propósitos. Estos atributos son objetos de tipo ndarray de NumPy, sin
embargo, lo que nos interesa no es el objeto como tal sino su contenido, el cual debe expresarse como un dato
de tipo float.
Para realizar la extracción de los valores correspondientes, se recurre al uso de la selección por medio de
corchetes “[ ]” y del método .item(), indicando la posición correspondiente al elemento de interés; asimismo,
emplearemos una Series de pandas para obtener una presentación un tanto más “agradable” del resultado. El
código a ejecutar es el siguiente:
73
Out[12]: Estadístico 55.4003
p-value 0.0000
dtype: float64
Asumiendo un nivel de significancia de 5 %, como p − value < α entonces se rechaza la hipótesis nula y se
concluye que los parámetros asociados a los instrumentos son estadísticamente distintos de cero (trabajando
con un α de 0.05). Se invita al lector a contrastar estos resultados con los presentados en Wooldridge (2010),
de modo que pueda comprobar que son exactamente iguales y que el proceso desarrollado es correcto.
Hasta este punto ya se ha llevado a cabo la regresión de la primera etapa y se ha verificado que los instru-
mentos empleados cumplen la condición de relevancia, por lo que es posible continuar con la segunda etapa
del método MC2E, la cual se aborda en la siguiente sección.
Segunda etapa
Durante la segunda etapa se regresa la variable dependiente de la ecuación original contra las variables ex-
plicativas exógenas y la variable estimada en la primera etapa (la cual reemplaza la variable explicativa endógena
que inicialmente causaba los inconvenientes) (Nótese que en esta etapa no se hace uso de los instrumentos).
Durante la primera etapa estimamos la variable explicativa endógena en función de las variables explicativas
exógenas y los instrumentos; sin embargo, aún desconocemos los valores de dicha variable, los cuales son los
que se emplearán en la segunda etapa. Para obtener tales valores, empleamos el método .predict() aplicado
a la instancia de resultados del modelo de la primera etapa (a la que hemos dado el conveniente nombre de
ResultadosPrimeraEtapa).
Los resultados generados por la aplicación del método .predict() serán asignados a una nueva variable
(que llamaremos predicted_educ) del DataFram data, de modo que podamos emplearlos con facilidad
posteriormente. El código a ejecutar es el siguiente:
Verificamos que la nueva variable haya sido efectivamente creada, recurriendo al método .head():
In [14]: [Link]()
Out[14]: inlf hours kidslt6 kidsge6 age educ wage repwage hushrs husage \
0 1 1610 1 0 32 12 3.3540 2.65 2708 34
1 1 1656 0 2 30 12 1.3889 2.65 2310 30
2 1 1980 1 3 35 12 4.5455 4.04 3072 40
3 1 456 0 3 34 12 1.0965 3.25 1920 53
4 1 1568 1 2 31 14 4.5918 3.60 2000 32
74
2 1.514138 225 12.771979
3 0.092123 36 11.767683
4 1.524272 49 13.914615
[5 rows x 23 columns]
Podemos observar que la última columna corresponde a la variable predicted_educ, la cual es la que
hemos creado previamente; los datos que contiene corresponden a los valores estimados durante la primera
etapa.
Ya contamos con todas las variables a emplear en la regresión de la segunda etapa. Ahora, al igual que
durante la primera etapa, agruparemos las regresoras correspondientes en un único objeto. Nuevamente, crea-
remos una lista Python “[...]” cuyo contenido corresponderá a las regresoras de la segunda etapa, es decir, las
variables explicativas exógenas del modelo original y la variable endógena estimada en la primera etapa; esta
lista se empleará como argumento de la función concat(...) de pandas (recuerde la asignación del valor de
1 al parámetro axis). Así:
Verificamos que el proceso se haya llevado a cabo exitosamente, recurriendo al método .head():
In [16]: [Link]()
Ahora, tal y como hemos señalado repetidamente, construimos el modelo tomando como variable depen-
diente la variable dependiente de la ecuación original (y1 , siguiendo la notación empleada a inicios del capítulo)
y como variables explicativas las variables explicativas exógenas de la ecuación original (las x, siguiendo la no-
tación inicial) y la variable estimada en la primera etapa (¡ŷ2 y no y2 !).
Dado que hemos agrupado las variables explicativas de la regresión de la segunda etapa en un único objeto
(es lo que hemos hecho en la última celda de código), basta con emplear dicho objeto como argumento en la
construcción del modelo (la cual se lleva a cabo con la ya conocida función OLS(...) de Statsmodels).
El código (que ya debe resultar familiar) para la construcción y estimación del modelo de la segunda etapa
es el siguiente:
75
Date: xxx, xx xxx xxxx Prob (F-statistic): 7.62e-05
Time: xx:xx:xx Log-Likelihood: -457.17
No. Observations: 428 AIC: 922.3
Df Residuals: 424 BIC: 938.6
Df Model: 3
Covariance Type: nonrobust
==================================================================================
coef std err t P>|t| [0.025 0.975]
----------------------------------------------------------------------------------
const 0.0481 0.420 0.115 0.909 -0.777 0.873
exper 0.0442 0.014 3.136 0.002 0.016 0.072
expersq -0.0009 0.000 -2.134 0.033 -0.002 -7.11e-05
predicted_educ 0.0614 0.033 1.863 0.063 -0.003 0.126
==============================================================================
Omnibus: 53.587 Durbin-Watson: 1.959
Prob(Omnibus): 0.000 Jarque-Bera (JB): 168.354
Skew: -0.551 Prob(JB): 2.77e-37
Kurtosis: 5.868 Cond. No. 4.41e+03
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
[2] The condition number is large, 4.41e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
Así concluye la aplicación del método MC2E en Python de forma manual. El lector debe tener en cuenta
que, aunque los parámetros estimados “manualmente” son correctos, los errores estándar no lo son (véase la
página 522 de Wooldridge (2010)), por lo que no es recomendable ejecutar las regresiones de primera y segun-
da etapa por cuenta propia. Para evitar este inconveniente, la mayoría de paquetes y programas estadísticos
ofrecen la posibilidad de ejecutar el método MC2E de forma automática, generando resultados correctos. El
lenguaje de programación Python no es la excepción, pues su librería linearmodels ofrece una función diseñada
específicamente para tal tarea, como veremos a continuación.
Proceso automático
Para evitar tener que hacer las regresiones de las dos etapas del método MC2E por nuestra propia cuenta,
la librería linearmodels de Python nos da la posibilidad de hacer uso de sus funciones para realizar “automática-
mente” dichas regresiones. La función específica empleada para tal propósito es IV2SLS(...) (acrónimo para
Instrumental Variables 2 Stages Least Squares), cuyos argumentos incluyen la variable dependiente (dependent), las
regresoras exógenas (exog), las regresoras endógenas (endog) y los instrumentos (instruments).
Para la estimación, como resulta familiar, se usa el método .fit(), y para el resumen de resultados el
atributo .summary.
El código correspondiente a nuestro ejemplo es:
76
ResultadosModelo2Etapas = [Link]()
print([Link])
Parameter Estimates
==============================================================================
Parameter Std. Err. T-stat P-value Lower CI Upper CI
------------------------------------------------------------------------------
const 0.0481 0.4278 0.1124 0.9105 -0.7903 0.8865
exper 0.0442 0.0155 2.8546 0.0043 0.0138 0.0745
expersq -0.0009 0.0004 -2.1001 0.0357 -0.0017 -5.997e-05
educ 0.0614 0.0332 1.8503 0.0643 -0.0036 0.1264
==============================================================================
Endogenous: educ
Instruments: motheduc, fatheduc
Robust Covariance (Heteroskedastic)
Debiased: False
El lector puede observar cómo la aplicación del método MC2E resulta tremendamente sencilla con el uso de
la función IV2SLS(...) de [Link] y que todas las tareas llevadas a cabo en el proceso manual han sido
ejecutadas perfectamente en tan solo ¡3 líneas de código!, sin tener que crear variables adicionales o emplear
múltiples grupos de las mismas, obteniendo, además, resultados completamente correctos.
77
OLS Regression Results
==============================================================================
Dep. Variable: lwage R-squared: 0.157
Model: OLS Adj. R-squared: 0.151
Method: Least Squares F-statistic: 27.56
Date: xxx, xx xxx xxxx Prob (F-statistic): 2.68e-16
Time: xx:xx:xx Log-Likelihood: -431.60
No. Observations: 428 AIC: 871.2
Df Residuals: 424 BIC: 887.4
Df Model: 3
Covariance Type: HC0
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
const -0.5220 0.201 -2.601 0.009 -0.915 -0.129
exper 0.0416 0.015 2.734 0.006 0.012 0.071
expersq -0.0008 0.000 -1.940 0.052 -0.002 8.28e-06
educ 0.1075 0.013 8.170 0.000 0.082 0.133
==============================================================================
Omnibus: 77.792 Durbin-Watson: 1.961
Prob(Omnibus): 0.000 Jarque-Bera (JB): 300.917
Skew: -0.753 Prob(JB): 4.54e-66
Kurtosis: 6.822 Cond. No. 2.21e+03
==============================================================================
Warnings:
[1] Standard Errors are heteroscedasticity robust (HC0)
[2] The condition number is large, 2.21e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
El lector puede observar que los resultados obtenidos por el método MCO difieren en cierta medida de los
correspondientes al método MC2E y que los errores estándar estimados por MC2E manualmente son diferentes
a los obtenidos automáticamente. Esto puede verificarlo al comparar las tablas de resumen de cada uno de los
modelos elaborados; sin embargo, tal ejercicio puede resultar tedioso, en cuanto cada una de estas tablas está
situada en una posición distinta dentro del notebook y comparar información entre estas puede resultar agota-
dor. Afortunadamente, la librería linearmodels brinda una herramienta para facilitar esta labor, como veremos
enseguida.
Comparación de modelos
Resulta muy práctico poder comparar los modelos estimados visualizando sus correspondientes resúme-
nes de instancias de resultados de forma simultánea; tal posibilidad puede materializarse gracias a la función
compare(...) que brinda la librería linearmodels en su módulo iv. Lo primero que debemos hacer para utilizar
la función compare(...) es contar con acceso a ésta, por lo que procedemos a importarla:
78
Habiendo ejecutado la celda de código anterior, habremos importado la función compare(...), por lo
que podremos hacer uso efectivo de la misma. Como argumento utilizaremos un diccionario Python “{...}”, en
el que emplearemos como claves los nombres que asignaremos a cada modelo y como valores las instancias de
resultados de los modelos correspondientes.
El usuario debe tener en cuenta que para que la función compare(...) de [Link] no genere un
error, todos los modelos a comparar deben corresponder a instancias de resultados del tipo generado por la
función IV2SLS(...). En este caso, como solo se utilizó la función IV2SLS(...) para el último modelo, no
es posible hacer una comparación con los demás. ¿Cómo puede solucionarse tal inconveniente? La respuesta es
muy simple en realidad: todos los modelos que hemos creado pueden construirse haciendo uso de la función
IV2SLS(...), tan solo hay que ser cuidadosos e incluir los argumentos correctos.
La función IV2SLS(...) tiene parámetros específicos para la variable dependiente, la(s) variable(s) exóge-
na(s), la(s) variable(s) endógena(s) y el(los) instrumento(s). Así, para una estimación común y corriente por el
método de MCO, basta con asignar None al parámetro endog y al parámetro instruments; para una estimación
por MC2E manual basta con emplear las regresoras de la segunda etapa en el parámetro exog y asignar None
a los parámetros endog e instruments; para una estimación por MC2E automática tan solo debe emplearse los
mismos argumentos que en el modelo elaborado previamente en la sección “Proceso automático”. El código
(se sugiere al lector examinarlo con detenimiento) para estimar los modelos correspondientes, empleando para
todos la función IV2SLS(...), es el siguiente:
Ahora, como todas las instancias de resultados han sido generadas por un modelo de la misma clase, po-
demos emplear la función compare(...) para comparar los resultados obtenidos y contrastar los modelos,
así:
Model Comparison
======================================================================
2 Etapas Automático 2 Etapas Manual MCO
----------------------------------------------------------------------
Dep. Variable lwage lwage lwage
Estimator IV-2SLS OLS OLS
No. Observations 428 428 428
Cov. Est. robust robust robust
R-squared 0.1357 0.0498 0.1568
Adj. R-squared 0.1296 0.0431 0.1509
F-statistic 18.611 17.111 82.671
P-value (F-stat) 0.0003 0.0007 0.0000
================== =========== =========== ===========
const 0.0481 0.0481 -0.5220
79
(0.1124) (0.1071) (-2.6010)
exper 0.0442 0.0442 0.0416
(2.8546) (2.7045) (2.7344)
expersq -0.0009 -0.0009 -0.0008
(-2.1001) (-1.9627) (-1.9402)
educ 0.0614 0.1075
(1.8503) (8.1697)
predicted_educ 0.0614
(1.7553)
==================== ============= ============= =============
Instruments motheduc
fatheduc
----------------------------------------------------------------------
La salida generada por la ejecución del código de la celda anterior nos permite visualizar de forma simultá-
nea los resultados de los diferentes modelos construidos, facilitando la comparación de los mismos. El lector
puede observar que en la parte inferior de la salida se indica T-stats reported in parentheses, lo que
señala que el número reportado entre paréntesis, debajo de cada uno de los coeficientes, corresponde al valor
t; sin embargo, en múltiples ocasiones lo que se reporta con correspondencia a un coeficiente no es el valor t
sino el error estándar. ¿Es posible reportar el error estándar en lugar del valor t? La respuesta es sí: la función
compare(...) posee un parámetro precision, con valor por defecto tstats, que permite especificar el estimador
de precisión (dentro de los disponibles) a incluir en la salida.
Para reportar errores estándar en lugar de valores t, basta con asignar “std_errors” al parámetro precision
de la función compare(...), como se evidencia a continuación:
Model Comparison
=====================================================================
2 Etapas Automático 2 Etapas Manual MCO
---------------------------------------------------------------------
Dep. Variable lwage lwage lwage
Estimator IV-2SLS OLS OLS
No. Observations 428 428 428
Cov. Est. robust robust robust
R-squared 0.1357 0.0498 0.1568
Adj. R-squared 0.1296 0.0431 0.1509
F-statistic 18.611 17.111 82.671
P-value (F-stat) 0.0003 0.0007 0.0000
================== ========== ========== ==========
const 0.0481 0.0481 -0.5220
(0.4278) (0.4492) (0.2007)
80
exper 0.0442 0.0442 0.0416
(0.0155) (0.0163) (0.0152)
expersq -0.0009 -0.0009 -0.0008
(0.0004) (0.0005) (0.0004)
educ 0.0614 0.1075
(0.0332) (0.0132)
predicted_educ 0.0614
(0.0350)
==================== ============ ============ ============
Instruments motheduc
fatheduc
---------------------------------------------------------------------
Nótese que los valores reportados entre paréntesis han cambiado y que en la sección inferior de la salida
aparece Std. Errors reported in parentheses en lugar del T-stats reported in parentheses
del caso previo. De este modo, lo que ahora se reporta en correspondencia a cada coeficiente es su error estándar
y no su valor t.
En este punto finaliza el tema de variables instrumentales y el método de Mínimos Cuadrados en 2 Etapas
en Python. El lector podrá observar que, hasta el momento, solo se ha trabajado con variables dependientes
continuas y que las regresiones llevadas a cabo en los distintos ejemplos son lineales; esto cambiará levemente
en el próximo capítulo, en donde se dará paso a una muy breve exposición de otro tipo de modelos.
81
Capítulo 6
En los capítulos previos se ha trabajado con modelos de regresión lineal empleando una variable continua
como variable dependiente; ahora, nos alejaremos un poco de los modelos con tales características, explorando
regresiones no lineales con variables dependientes discretas.
En los primeros modelos que abordaremos, la variable dependiente es una variable categórica, con solo dos
categorías, codificada con los valores 0 y 1. Dicha variable solo toma estos dos valores, los cuales no tienen un
significado por sí mismos (a diferencia del caso de variables continuas), en cuanto su función es de codificadores:
toma el valor de 1 cuando la observación pertenece a un grupo específico y de 0 cuando no pertenece a tal grupo.
En palabras un tanto más simples: la variable toma el valor de 1 cuando se cumple una condición (pertenencia
a una clasificación específica) y de 0 cuando no se cumple dicha condición.
Es posible que el lector no tenga completa claridad sobre las características de las variables que emplea-
remos y que la descripción anterior no le haya resultado del todo satisfactoria; en tal caso, puede que algunos
ejemplos sean de ayuda:
Si contamos con un grupo formado por personas de ambos sexos, en dicho grupo habrá hombres y ha-
brá mujeres; si empleamos una variable ’femenino’ para realizar una clasificación por sexo, esta variable
tomará el valor de 1 cuando la información corresponda a una mujer y de 0 cuando corresponda a un
hombre. ’femenino’ es una variable binaria, pues toma un valor cuando una condición se cumple (1 cuan-
do la condición de sexo femenino se cumple) y toma otro valor cuando la condición no se cumple (0
cuando la condición de sexo femenino no se cumple).
Si tenemos un grupo formado por hombres, es posible que algunos de estos hombres estén casados
mientras que otros no, así, podemos emplear una variable binaria ’casado’ para clasificarlos en el grupo
correspondiente. La variable ’casado’ toma el valor de 1 cuando corresponde a un hombre casado y toma
el valor de 0 cuando corresponde a un hombre no casado.
En un grupo de estudiantes universitarios algunos de estos habrán perdido asignaturas mientras que
otros no. De este modo, podemos emplear una variable binaria que tome el valor de 1 cuando corres-
ponda a un estudiante que ha perdido asignaturas y de 0 cuando corresponda a un estudiante que no ha
perdido asignaturas.
El lector puede observar que la variable binaria toma el valor de 1 cuando se cumple la condición específica
y de 0 en caso contrario; tal asignación de valores no tiene sentido numérico y no afecta la regresión. Tranquila-
mente se puede asignar el valor de 0 cuando la condición se cumple y de 1 en caso contrario; sin embargo, por
convención, y mayor facilidad en la interpretación, se emplea 1 para el cumplimiento de la condición y 0 para
el no cumplimiento de la misma.
82
Hasta este punto el lector tiene una idea general de lo que es una variable binaria, sin embargo, aún no se le
ha informado cómo se empleará dicha variable en la regresión. Es muy importante señalar que en los modelos
que abordaremos a continuación, lo que se estima en realidad es la probabilidad de que la variable dependiente
tome el valor de 1 (al lector interesado en conocer el porqué, se recomienda consultar el Capítulo 15, Modelos
de regresión de respuesta cualitativa, de Gujarati y Porter (2010)); dicha probabilidad está dada por los parámetros
(β) y las variables explicativas (x).
En términos un tanto más técnicos:
Prob[yi = 1| xi ] = Pi = F ( xiT β)
en donde xiT corresponde al vector transpuesto de variables explicativas y β al vector de parámetros. Así,
los coeficientes que se obtienen durante la estimación se asocian al impacto que tiene determinada variable
explicativa sobre la probabilidad de que la variable dependiente tome el valor de 1; en otras palabras, los coefi-
cientes se relacionan (no necesariamente cuantifican de forma directa, debido a la forma funcional adoptada,
como veremos más adelante) a la variación en la probabilidad de que la variable dependiente tome el valor de
1.
Algunos de los modelos empleados para trabajar con variables de respuesta binaria son: 1) Modelo Lineal
de Probabilidad, 2) Modelo Logit y 3) Modelo Probit. Es posible construir y estimar fácilmente cada uno de
estos modelos en Python con la ayuda de librerías especializadas. Tales labores son las que llevaremos a cabo a
continuación, no sin antes realizar una muy mínima exposición de las ideas detrás de cada uno de los modelos
tratados.
El dataset que utilizaremos para desarrollar los distintos modelos será el empleado en el Ejemplo 15.7,
Fumar o no fumar, de Gujarati y Porter (2010). Este conjunto de información se ha obtenido del sitio web eco-
[Link] de SHAZAM Analytics Ltd.
83
import [Link] as sm
import [Link] as sms
import [Link] as plt
El anterior bloque de código permite cargar las herramientas necesarias para llevar a cabo las acciones
requeridas en el análisis econométrico planteado; en caso de que el lector no tenga claridad al respecto, se le
pide consultar la sección “Preparación del entorno” del Capítulo 2, en donde se brinda una explicación sobre
tal cuestión.
84
Podemos observar que el proceso de importación ha sido exitoso y la información del dataset ha sido alma-
cenada correctamente en un DataFrame de pandas. Ahora, podemos continuar tranquilamente con el desarrollo
de nuestro análisis, en cuanto hemos concluido satisfactoriamente el primer paso.
85
El lector que haya seguido esta obra hasta este punto habrá notado que en los capítulos previos hemos
almacenado diversas variables en objetos específicos. Esto lo hemos hecho para agrupar dichas variables de
acuerdo a su rol particular en la regresión, tratando de simplificar el código a utilizar posteriormente. Ahora,
para los distintos modelos que abordaremos en este capítulo se utilizarán las mismas variables, agrupadas del
mismo modo para todos. Procederemos a realizar la agrupación por medio de los corchetes “[ ]” y las listas
Python “[...]”, tal y como hemos hecho en capítulos previos:
In [6]: Y = data["Fumador"]
X = data[["Edad", "Escolaridad", "Ingreso", "Pcigs79"]]
El primer modelo a examinar es el Modelo Lineal de Probabilidad, el cual es, básicamente, un modelo de
regresión lineal con variable dependiente binaria y estimación por el método de Mínimos Cuadrados Ordinarios.
La regresión lineal y el método MCO se trataron en el Capítulo 3 de esta obra, por lo que a quien no tenga claridad
al respecto se le sugiere consultar dicho capítulo.
Para la construcción del modelo y estimación de los coeficientes se emplean (como lo hemos hecho en oca-
siones previas) la función OLS(...) de Statsmodels y el método .fit(); para la visualización de los resultados,
la función print(...) y el método .summary(). El código a ejecutar es el siguiente:
86
OLS Regression Results
==============================================================================
Dep. Variable: Fumador R-squared: 0.039
Model: OLS Adj. R-squared: 0.036
Method: Least Squares F-statistic: 12.01
Date: xxx, xx xxx xxxx Prob (F-statistic): 1.43e-09
Time: xx:xx:xx Log-Likelihood: -809.19
No. Observations: 1196 AIC: 1628.
Df Residuals: 1191 BIC: 1654.
Df Model: 4
Covariance Type: nonrobust
===============================================================================
coef std err t P>|t| [0.025 0.975]
-------------------------------------------------------------------------------
const 1.1231 0.188 5.963 0.000 0.754 1.493
Edad -0.0047 0.001 -5.701 0.000 -0.006 -0.003
Escolaridad -0.0206 0.005 -4.465 0.000 -0.030 -0.012
Ingreso 1.026e-06 1.63e-06 0.629 0.530 -2.18e-06 4.23e-06
Pcigs79 -0.0051 0.003 -1.799 0.072 -0.011 0.000
==============================================================================
Omnibus: 37.223 Durbin-Watson: 1.944
Prob(Omnibus): 0.000 Jarque-Bera (JB): 173.523
Skew: 0.449 Prob(JB): 2.09e-38
Kurtosis: 1.364 Cond. No. 2.91e+05
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly
specified.
[2] The condition number is large, 2.91e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
Uno de los aspectos más importantes a considerar en los modelos de variable dependiente binaria es el
efecto marginal de cada una de las variables explicativas; tal efecto corresponde a la variación en la probabili-
dad (de que la variable dependiente tome el valor de 1) respecto a la variación en una variable explicativa. En
∂F ( x T β)
términos un poco más técnicos, el efecto marginal de una variable x j corresponde al valor de ∂xiji , en donde
x j es una variable que hace parte del vector x.
Así, para el Modelo Lineal de Probabilidad, los efectos marginales no son más que los coeficientes, pues:
∂F ( xiT β)
= βj
∂x ji
Por lo tanto, para conocer los efectos marginales del Modelo Lineal de Probabilidad en Python, basta con
conocer los coeficientes asociados a las variables. Estos coeficientes corresponden a uno de los atributos (que
pueden consultarse con la función dir(...)) de la instancia de resultados, y se identifica como .params.
87
Podemos asignar la aplicación de este atributo a un objeto específico y visualizar su contenido utilizando la
función print(...). Así:
Así, observamos, por ejemplo, que un incremento de un año en la edad está asociado, en promedio, a una
disminución de 0.47 puntos porcentuales en la probabilidad de ser fumador, y que un aumento de 1 año en la
escolaridad está asociado, en promedio, a una disminución de 2 puntos porcentuales en la probabilidad de ser
fumador.
En realidad, para nuestros propósitos, no hay más aspectos a tratar con detenimiento en el Modelo Lineal
de Probabilidad; basta con conocer su proceso de construcción y estimación y la obtención de los efectos mar-
ginales. Ahora, tal y como se señaló a inicios del capítulo, nos alejaremos de los modelos de regresión lineal y
exploraremos otras posibilidades: la primera, el Modelo Logit.
Modelo Logit
El Modelo Lineal de Probabilidad, aunque muy simple, no tiende a ser ampliamente usado en la práctica;
la razón para tal situación se encuentra, principalmente, en el hecho de que este modelo puede llevar a la
obtención de probabilidades carentes de sentido (por ejemplo, mayores a 1 o que violan el axioma de que la
probabilidad no puede ser negativa).
Aunque no se mencionó en la sección anterior, el MLP puede resultar muy problemático, en cuanto puede
generar probabilidades menores a 0 o mayores a 1, lo que le resta atractivo y conveniencia. Para superar las
deficiencias del Modelo Lineal de Probabilidad se recurre a modelos de regresión no lineales, de los cuales los
más populares son el Logit y el Probit.
En el Modelo Logit, la probabilidad de que la variable dependiente tome el valor de 1 se expresa como una
función de la siguiente forma:
1
Pi = F xiT β =
1 + e−( β0 + β1 x1i +...+ β k xki )
Alternativamente, tomando Zi = β 0 + β 1 x1i + ... + β k xki , la función puede expresarse de la siguiente
manera:
1
Pi = F ( Zi ) =
1 + e−Zi
e Zi 1
Otro modo de escribir la expresión anterior es Pi = 1+e Zi
, por lo que 1 − Pi = 1+e Zi
. Así, se tiene que:
Pi
= e Zi
1 − Pi
88
Aplicando el logaritmo a la anterior expresión se obtiene:
Pi
L = ln = Zi = β 0 + β 1 x1i + ... + β k xki
1 − Pi
L se denomina logit y corresponde al logaritmo de la razón de las probabilidades. Así, “mientras el MLP
supone que Pi está linealmente relacionado con xi , el modelo logit supone que el logaritmo de la razón de
probabilidades está relacionado linealmente con xi ” (Gujarati y Porter, 2010, p. 555).
La estimación de los parámetros (β) en el Modelo Logit es llevada a cabo por medio del método de Máxima
Verosimilitud (Maximum Likelihood) o MV. Al lector interesado en conocer con mayor detalle el Modelo Logit
y su proceso de estimación, se sugiere consultar el Capítulo 15, Modelos de regresión de respuesta cualitativa, de
Gujarati y Porter (2010).
Uno de los puntos a considerar en el Modelo Logit es que la interpretación de los coeficientes no resulta de
mayor utilidad, en cuanto estos cuantifican el impacto que tiene el cambio en una variable explicativa sobre el
logaritmo de la razón de probabilidades, algo que, a primera vista, no tiene gran valor ilustrativo. Debido a esto, se
recurre a los efectos marginales para cuantificar el impacto de las variables explicativas.
Como resulta evidente, el efecto marginal de una variable explicativa en un Modelo Logit, a diferencia del
caso del Modelo Lineal de Probabilidad, no corresponde al coeficiente asociado a dicha variable, sino que está
dado por:
∂F xiT β
e Zi
= βj 2
∂x ji (1 + e Zi )
Así, se observa que el efecto marginal de una variable explicativa depende del conjunto de variables expli-
cativas y no está dado directamente por el coeficiente asociado a ésta. Asimismo, es importante notar que en la
expresión anterior se ha hallado el efecto marginal de la variable x j para la observación i, pero debe tenerse en
consideración que el dataset contiene información de múltiples observaciones y los valores de las variables son
diferentes entre dichas observaciones, por lo que no existe un único valor para el efecto marginal de la variable
x j sino n (número de observaciones) valores para éste.
Dado que no resulta muy práctico conocer la magnitud del efecto marginal de una variable para cada una de
las observaciones, es necesario contar con una medida “general” del efecto marginal de determinada variable,
pero ¿cómo se obtiene dicha medida “general”? Una respuesta posible a este interrogante es: existen varios
“caminos”. A continuación, se presentan los más populares.
El primer camino (conocido como efectos marginales en la media) para obtener una medida “general” del
efecto marginal de una variable es utilizar la observación promedio para calcularlo; es decir, emplear el promedio
de cada variable como valor de la variable correspondiente dentro de Z.
Definiendo el vector que contiene los valores promedio de las variables explicativas como x̄, y teniendo en
cuenta que, consecuentemente, x̄ T β = Z̄, el efecto marginal de la variable x j está dado por:
∂P e Z̄
= βj 2
∂x j 1 + e Z̄
Así, el efecto marginal “general” de determinada variable corresponde al efecto marginal de dicha variable
calculado con los valores correspondientes a la observación promedio.
El otro camino (conocido como efectos marginales promedio) para obtener una medida “general” del efecto
marginal de una variable es obtener el promedio de los efectos marginales individuales correspondientes a dicha
variable; es decir, calcular el efecto marginal de la variable x j para cada una de las observaciones i (en total se
tienen n observaciones), y hallar el promedio de dichos efectos marginales. Concretamente, el efecto marginal
de la variable x j viene dado por:
89
∂F xiT β e Zi
∑in=1
∂P ∑in=1 ∂x ji 1+e Zi
2
= = βj
∂x j n n
Así, el efecto marginal “general” de determinada variable corresponde al promedio de los efectos marginales
individuales para dicha variable.
Hasta este punto se ha dado una descripción, muy breve, de la idea general sobre la que trata el Modelo
Logit; ahora, es momento de aplicar a la práctica los conocimientos adquiridos, llevando a cabo la construcción
y estimación de un modelo logit en Python.
90
Escolaridad -0.0910 0.021 -4.402 0.000 -0.131 -0.050
Ingreso 4.72e-06 7.17e-06 0.658 0.510 -9.33e-06 1.88e-05
Pcigs79 -0.0223 0.012 -1.789 0.074 -0.047 0.002
===============================================================================
El lector debe recordar que la interpretación directa de los coeficientes en el Modelo Logit no resulta par-
ticularmente útil a primera vista, en cuanto estos coeficientes cuantifican el impacto de determinada variable
sobre el logaritmo de la razón de probabilidades y no directamente sobre la probabilidad de que la variable depen-
diente tome el valor de 1.
Para conocer el impacto de una variable específica en un modelo logit, tal y como se señaló a ini-
cio de este capítulo, se recurre a los efectos marginales. Estos efectos se obtienen por medio del método
.get_margeff(), aplicado a la instancia de resultados del modelo. El lector debe tener en cuenta que pa-
ra obtener una medida “general” del efecto marginal de una variable específica se señaló la existencia de dos
caminos (que no son los únicos): 1) obtener el efecto marginal usando la observación promedio y 2) obtener el
promedio de los efectos marginales individuales.
Para obtener los efectos marginales utilizando el primer “camino”, basta con asignar el valor ’mean’ al
parámetro at del método .get_margeff(); para obtenerlos por el segundo “camino”, basta con asignarle el
valor de ’overall’. El parámetro at del método .get_margeff() corresponde a ’overall’ por defecto, por lo
que los efectos marginales promedio son los que se obtienen en caso de que no se especifique el valor de at.
Los efectos marginales promedio (segundo camino) son:
In [10]: LogitMargEff = ResultadosLogit.get_margeff()
print([Link]())
Logit Marginal Effects
=====================================
Dep. Variable: Fumador
Method: dydx
At: overall
===============================================================================
dy/dx std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
Edad -0.0047 0.001 -5.864 0.000 -0.006 -0.003
Escolaridad -0.0206 0.005 -4.539 0.000 -0.030 -0.012
Ingreso 1.069e-06 1.62e-06 0.659 0.510 -2.11e-06 4.25e-06
Pcigs79 -0.0051 0.003 -1.798 0.072 -0.011 0.000
===============================================================================
Esta tabla, además del valor específico de los efectos marginales, también presenta los intervalos de con-
fianza correspondientes e indica qué método es empleado para calcularlos (At:). Así, observamos, por ejemplo,
que un incremento de un año en la edad está asociado, en promedio, a una disminución de 0.47 puntos por-
centuales en la probabilidad de ser fumador y que un aumento de 1 año en la escolaridad está asociado, en
promedio, a una disminución de 2 puntos porcentuales en la probabilidad de ser fumador.
Para los efectos marginales en la media (primer camino), tan solo se necesita asignar ’mean’ al parámetro
at. El código es el siguiente:
In [11]: LogitMargEff = ResultadosLogit.get_margeff(at = "mean")
print([Link]())
91
Logit Marginal Effects
=====================================
Dep. Variable: Fumador
Method: dydx
At: mean
===============================================================================
dy/dx std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
Edad -0.0049 0.001 -5.598 0.000 -0.007 -0.003
Escolaridad -0.0213 0.005 -4.411 0.000 -0.031 -0.012
Ingreso 1.107e-06 1.68e-06 0.658 0.510 -2.19e-06 4.4e-06
Pcigs79 -0.0052 0.003 -1.790 0.073 -0.011 0.000
===============================================================================
Así, observamos, por ejemplo, que un incremento de un año en la edad está asociado, en promedio, a una
disminución de 0.49 puntos porcentuales en la probabilidad de ser fumador y que un aumento de 1 año en la
escolaridad está asociado, en promedio, a una disminución de 2.13 puntos porcentuales en la probabilidad de
ser fumador.
Modelo Probit
Al igual que el Modelo Logit, una alternativa popular al MLP es el Modelo Probit. En este modelo, la pro-
babilidad de que la variable dependiente tome el valor de 1 viene dada por lo siguiente:
Z xT
i
xiT β =Φ xiT β
Pi = F = φ(z)dz
−∞
igual que con el Modelo Logit, las probabilidades obtenidas no son inferiores a 0 o superiores a 1 y el método de
estimación de los parámetros es el de Máxima Verosimilitud (MV). Al lector interesado en conocer con mayor
detalle el Modelo Probit y su proceso de estimación, se le recomienda consultar el Capítulo 15, Modelos de
regresión de respuesta cualitativa, de Gujarati y Porter (2010).
Como en el caso del Modelo Logit, la interpretación directa de los coeficientes asociados a las variables
explicativas no resulta de gran utilidad y los efectos marginales correspondientes difieren de estos. Para el
Modelo Probit, el efecto marginal de la variable x j viene dado por:
∂F xiT β
= φ xiT β β j
x ji
en donde φ(.) es la Función de Densidad de Probabilidad de la Distribución Normal Estándar. Como puede
observarse, los efectos marginales en el Modelo Probit, al igual que en el Modelo Logit, no están dados directa-
mente por el coeficiente asociado a la variable, sino que dependen de los valores de las variables explicativas.
Así, en un dataset con n observaciones se obtiene n valores para el efecto marginal de la misma variable; para
obtener una medida “general”, pueden seguirse los dos caminos previamente señalados en el caso del Modelo
Logit.
En el primer camino se emplean los valores promedio de las variables con el propósito de hallar un único
efecto marginal para determinada variable. Definiendo, nuevamente, el vector que contiene los valores prome-
dio de las variables explicativas como x̄, el efecto marginal de la variable x j está dado por:
92
∂P
= φ x̄iT β
∂x j
Así, el efecto marginal ’general’ de determinada variable corresponde al efecto marginal de dicha variable
calculado con los valores correspondientes a la observación promedio.
En el segundo camino, se calculan los n efectos marginales correspondientes a la misma variable y simple-
mente se obtiene su promedio. Concretamente:
∑in=1 φ xiT β
∂P
= βj
∂x j n
Así, el efecto marginal “general” de determinada variable corresponde al promedio de los efectos marginales
individuales para dicha variable.
Hasta este punto se ha dado una descripción muy breve de la idea general sobre la que trata el Modelo
Probit; ahora, es momento de aplicar a la práctica los conocimientos adquiridos, llevando a cabo la construcción
y estimación de un modelo probit en Python.
93
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
const 1.7019 0.511 3.333 0.001 0.701 2.703
Edad -0.0130 0.002 -5.655 0.000 -0.017 -0.008
Escolaridad -0.0562 0.013 -4.450 0.000 -0.081 -0.031
Ingreso 2.72e-06 4.4e-06 0.619 0.536 -5.9e-06 1.13e-05
Pcigs79 -0.0138 0.008 -1.792 0.073 -0.029 0.001
===============================================================================
El lector debe recordar que la interpretación directa de los coeficientes en el Modelo Probit no resulta
particularmente útil a primera vista y que para conocer el impacto de una variable específica en un modelo
probit, tal y como se señaló a inicio de este capítulo, se recurre a los efectos marginales.
Los efectos marginales se obtienen por medio del método .get_margeff(), aplicado a la instancia de re-
sultados del modelo. El lector debe tener en cuenta que para obtener una medida “general” del efecto marginal
de una variable específica se señaló la existencia de dos caminos (que no son los únicos): 1) obtener el efecto
marginal usando la observación promedio y 2) obtener el promedio de los efectos marginales individuales.
Para obtener los efectos marginales utilizando el primer “camino” basta con asignar el valor ’mean’ al pará-
metro at del método .get_margeff(); para obtenerlos por el segundo “camino”, basta con asignarle el valor
de ’overall’. El parámetro at del método .get_margeff() corresponde a ’overall’ por defecto, por lo que los
efectos marginales promedio son los que se obtienen en caso de que no se especifique el valor de at.
Los efectos marginales promedio (segundo camino) son:
La información que se presenta en esta tabla es la misma que se obtiene en el caso del modelo Logit. Así, ob-
servamos, por ejemplo, que un incremento de un año en la edad está asociado, en promedio, a una disminución
de 0.48 puntos porcentuales en la probabilidad de ser fumador y que un aumento de 1 año en la escolaridad
está asociado, en promedio, a una disminución de 2 puntos porcentuales en la probabilidad de ser fumador.
Para los efectos marginales en la media (primer camino) solo es necesario modificar el valor del parámetro
at; el código es el siguiente:
94
In [14]: ProbitMargEff = ResultadosProbit.get_margeff(at = "mean")
print([Link]())
Empleando este método, observamos, por ejemplo, que un incremento de un año en la edad está asociado,
en promedio, a una disminución de 0.49 puntos porcentuales en la probabilidad de ser fumador y que un au-
mento de 1 año en la escolaridad está asociado, en promedio, a una disminución de 2.13 puntos porcentuales
en la probabilidad de ser fumador.
Comparación de modelos
Resulta útil poder comparar los resultados de los distintos modelos tratados teniendo a disposición su
información en un único espacio; es decir, contando con una tabla de resumen en la que se presente simultá-
neamente la información básica de la instancia de resultados correspondiente a cada modelo.
Es posible comparar resultados de modelos creados con funciones de Statsmodels empleando la función
summary_col(...) del módulo iolib.summary2 de Statsmodels. Esta función permite resumir múltiples
instancias de resultados presentándolas una al lado de la otra; algunos de sus argumentos son la lista de ins-
tancias de resultados a comparar (results), la lista de nombres a asignar al resumen de cada instancia (model_na-
mes) y el parámetro stars (de tipo booleano), que indica si se quiere o no visualizar la significancia por medio de
asteriscos (*).
El primer paso a seguir para poder comparar las instancias de resultados de los modelos es importar la fun-
ción que permite llevar a cabo tal tarea. La línea de importación de la función summary_col(...) de Statsmo-
[Link].summary2 es la siguiente:
Ahora, teniendo a disposición la función requerida, basta con especificar los argumentos correspondientes
y ejecutar el código para obtener la tabla de resumen (que denominaremos MiTabla):
95
=================================================
Modelo MLP Modelo Logit Modelo Probit
-------------------------------------------------
const 1.1231*** 2.7451*** 1.7019***
(0.1884) (0.8292) (0.5106)
Edad -0.0047*** -0.0209*** -0.0130***
(0.0008) (0.0037) (0.0023)
Escolaridad -0.0206*** -0.0910*** -0.0562***
(0.0046) (0.0207) (0.0126)
Ingreso 0.0000 0.0000 0.0000
(0.0000) (0.0000) (0.0000)
Pcigs79 -0.0051* -0.0223* -0.0138*
(0.0029) (0.0125) (0.0077)
=================================================
Standard errors in parentheses.
* p<.1, ** p<.05, ***p<.01
El lector debe recordar que la interpretación directa de los coeficientes de los modelos Logit y Probit no
resulta particularmente útil y que no es posible comparar su magnitud con la de los correspondientes al Modelo
Lineal de Probabilidad, en tanto cuantifican relaciones distintas. Debido a esto, se ha recurrido a los efectos
marginales para poder obtener una medida comparable, entre modelos, del efecto de una variable explicativa
sobre la probabilidad de que la variable dependiente tome el valor de 1. La tabla de resumen que hemos apenas
creado presenta los coeficientes de cada modelo y no los efectos marginales, ¿esto quiere decir que no es posible
comparar la información presentada?
Aunque no es posible comparar los coeficientes de los modelos MLP, Logit y Probit directamente y, por
tanto, utilizamos los efectos marginales, el lector debe tener en cuenta que existe una relación aproximada
entre las magnitudes de dichos coeficientes.
Siguiendo a Katchova (2013), los coeficientes de los modelos MLP, Logit y Probit se caracterizan por las
siguientes relaciones:
96
El objetivo que pretendemos es obtener una tabla de resumen en la que se presenten los efectos marginales
correspondientes al Modelo Lineal de Probabilidad, al Modelo Logit y al Modelo Probit, por lo que los datos
indispensables que necesitamos son los valores de tales efectos marginales. Lo primero que haremos es crear
DataFrames de pandas para almacenar los datos correspondientes.
La información de los efectos marginales del Modelo Logit y del Modelo Probit se ha almacenado con an-
terioridad en objetos específicos; ahora, extraeremos únicamente las magnitudes de dichos efectos marginales
(con el atributo .margeff, aplicado al objeto que contiene la información de los efectos marginales respecti-
vos) y las emplearemos como valores del parámetro data en la función DataFrame(...) de pandas. Asimismo,
para el parámetro index haremos uso de una lista Python “[...]” que contenga los nombres de las variables co-
rrespondientes y asignaremos al parámetro columns el nombre del respectivo modelo. El procedimiento es el
mismo para los modelos Logit y Probit, como se evidencia a continuación:
Ya contamos con los datos de los efectos marginales de los modelos Logit y Probit, pero, ¿qué hay de los
efectos marginales del MLP? El lector debe recordar que los efectos marginales del Modelo Lineal de Probabi-
lidad corresponden directamente a los coeficientes, por lo que, a diferencia de los modelos Logit y Probit, su
obtención se logra con la aplicación del atributo .params. El resultado de la aplicación de este atributo difiere
del que se obtiene con la aplicación del método .get_margeff() de los modelos Logit y Probit, por lo que
el procedimiento para extraer las magnitudes de los efectos marginales (coeficientes) será un tanto distinto al
que hemos desarrollado previamente.
Lo primero que el usuario debe notar es que el atributo .params genera un resultado en el que se incluye
el coeficiente correspondiente al término del intercepto, sin embargo, para los modelos Logit y Probit, como re-
sulta lógico, no se reporta un efecto marginal correspondiente a éste. Por lo tanto, con el propósito de comparar
los efectos marginales de las mismas variables, se ignorará el término del intercepto (β 0 ).
Para realizar la “extracción” de las magnitudes de todos los coeficientes exceptuando el término del inter-
cepto se hace uso de los corchetes “[ ]”, aplicados al objeto que contiene los parámetros del MLP (previamente
lo almacenamos como una Series de pandas). Dentro de los corchetes se empleará la secuencia 1:5; esto,
con el propósito de ignorar el término del intercepto. Adicionalmente, el nombre del modelo se asignará al
parámetro columns.
Nota: Al lector que lo desconozca, se le informa que una Series de pandas indiza (por defecto) sus ob-
servaciones con números enteros desde 0 y que la selección por corchetes “[ ]” con enteros es incluyente
en el primer valor y excluyente en el último, por lo que el empleo de una secuencia 1:5 dentro de cor-
chetes “[ ]” de selección significa “ignorar la primera observación (que está identificada con 0) e incluir
todas las demás hasta la indizada con el número 4 (no incluye la que tiene 5 como índice)”.
Ya se cuenta con los DataFrame que contienen las magnitudes de los efectos marginales de los modelos
Logit, Probit y MLP. Ahora, es momento de agrupar estos DataFrames en uno solo; para lograrlo, se emplea
la función concat(...) de pandas; en sus argumentos se incluye una lista Python “[...]” que contiene los
DataFrame a concatenar y el parámetro axis con valor de 1, de modo que la concatenación se logre ubicando
un DataFrame al lado de otro:
97
In [19]: [Link]([EfectosMLP, EfectosLogit, EfectosProbit], axis = 1)
Ahora se cuenta con una única tabla de resumen en la que se presentan los efectos marginales de los distin-
tos modelos: magnitudes que sí son directamente comparables. Podemos observar que los efectos estimados
no difieren significativamente: por ejemplo, el efecto marginal de la variable Escolaridad es de -2.06 p.p en el
MLP, de -2.13 p.p en el modelo Logit y de -2.13 p.p en el modelo Probit; para la variable Edad el efecto marginal
va desde -0.4726 p.p (MLP) hasta -0.492 p.p (Probit).
En este punto cierra este capítulo, en el cual hemos tratado los modelos con variable dependiente binaria;
ahora, extenderemos un tanto el alcance de los modelos de variable dependiente discreta al incluir regresadas
que no toman únicamente dos valores, es decir, variables no binarias.
98
Capítulo 7
En el capítulo previo se trabajó con modelos en los que la variable dependiente era discreta y, además, solo
podía tomar dos valores; en este capítulo, se dará continuidad al tema de la variable dependiente discreta, pero
se eliminará la restricción de que ésta solo puede tomar dos valores, es decir, abordaremos el tema de variables
discretas no binarias.
Las variables dependientes con las que trabajaremos serán categóricas, pero no estarán limitadas a solamen-
te dos categorías, por lo que no tomarán únicamente dos valores (0 y 1) sino tantos valores como alternativas
tenga dicha variable.
La breve exposición de los fundamentos teóricos del modelo a desarrollar en este capítulo está basada, en
su totalidad, en Katchova (2013, 2), por lo que, en caso de ser necesario, se invita al lector a consultar dicha
fuente para tener mayor claridad al respecto.
Las variables dependientes a emplear son discretas no binarias, por lo que toman valores de un conjunto
finito compuesto por más de dos opciones. Es posible que algunos ejemplos tengan cierto valor para ilustrar la
naturaleza de las mencionadas variables:
Cuando una persona desea consumir una cerveza, tiene a su disposición un conjunto finito y bien defi-
nido de marcas de cerveza que puede adquirir. Estas marcas pueden codificarse por medio de números,
asignando un número distinto a cada una, por lo que la variable resultante corresponde a una variable dis-
creta no binaria (asumiendo que existen más de dos marcas de cerveza a disposición), en cuanto existen
múltiples, pero limitadas, posibilidades.
Cuando alguien asiste a una sala de cine, tiene a su disposición distintos “combos” de comida. Existe un
conjunto bien definido de “combos” y su cantidad es limitada; de este modo, se puede asignar un número
distinto a cada combo, con el propósito de identificarlo. Así, la variable resultante de la codificación de
los “combos” es una variable discreta no binaria (asumiendo que existen más de dos combos).
A una variable discreta con m categorías (o alternativas) se aplica un proceso de codificación empleando
m números; estos números no tienen significado por sí mismos, en cuanto su magnitud no es interpretable y
su función es exclusivamente de codificadores. Asimismo, debe notarse que el orden en el que se asignan los
números es libre, en cuanto estos carecen de significado propio y no tienen naturaleza ordinal, por lo que no
establecen un orden definido para las alternativas de las variables.
La variable tiene múltiples categorías, pero es importante señalar que a determinado individuo (entidad)
solo corresponde una única alternativa; es decir, cada individuo o entidad, de las múltiples alternativas a dis-
posición, puede elegir únicamente una. Retomando los ejemplos anteriores, para el caso de las cervezas, el
individuo solo puede elegir una marca para la cerveza que consumirá en este momento y, para el caso de los
combos de comida, la persona solo debe elegir un único combo.
99
Para este tipo de modelos, la información se presenta en dos formatos específicos: wide y long. En el formato
wide, la información para cada individuo (entidad) se presenta en una única fila, por lo tanto y = j ( j =
1, 2, 3, ..., m); en el formato long, la información para cada individuo (entidad) se presenta en m filas, cada una
de las cuales corresponde a una de las alternativas de la variable dependiente, la cual toma el valor de 1 una
única vez y de 0 m − 1 veces para cada individuo (entidad).
Recurriendo a una descripción un tanto más ilustrativa, a continuación, se presenta la misma información,
en formato wide y en formato long. Retomando el ejemplo de las cervezas, tomando en consideración solo 5
marcas (codificadas de la siguiente manera: Aguila : 1, Poker : 2, [Link] : 3, Corona : 4, Heineken :
5), la siguiente información se presenta en formato wide:
Así, el lector puede apreciar que para el formato wide y = j, mientras que para el formato long:
100
(
1, si y = j
yj =
0, si y 6= j
Algo importante que debe mencionarse es que en los modelos de variable dependiente discreta no binaria,
al igual que en los modelos de variable dependiente discreta binaria, lo que interesa no es el valor de la variable
dependiente como tal sino la probabilidad de que esta tome un valor en específico. Así, en el Modelo Logit,
puesto que la variable dependiente era binaria, solo existían dos posibilidades, por lo que se estimaba la proba-
bilidad de que la variable dependiente tomará el valor de 1, y la probabilidad de que esta tomará el valor de 0
simplemente correspondía a uno menos la probabilidad de que tomará el valor de 1. Sin embargo, en el modelo
Logit Multinomial (que es el que trabajaremos en adelante) no existen únicamente dos alternativas, entonces,
¿qué probabilidad es la que se estima?
El lector debe saber que en el modelo Logit Multinomial lo que se busca es estimar la probabilidad de que se
elija la alternativa j dentro de las m alternativas que tiene la variable dependiente. Para los modelos de variable
dependiente discreta no binaria, la probabilidad de que el individuo i elija la alternativa j está dada por:
101
en donde Pji es la probabilidad de que el individuo (entidad) i elija la alternativa j, wi es el vector de
variables alternative-invariant para el individuo (entidad) i y γ j es el conjunto de coeficientes correspondientes
a la alternativa j. Para lograr la estimación de la probabilidad, uno de los conjuntos de coeficientes se normaliza
a 0, de modo que la interpretación de los demás se realiza con referencia a la categoría base (a la que corresponde
el conjunto de coeficientes 0).
Un punto a destacar es que la magnitud de los coeficientes del Modelo Logit Multinomial, al igual que de
los Logit y Probit, no debe interpretarse directamente; estos coeficientes tan solo indican si aumenta (cuando
el coeficiente es positivo) o disminuye (cuando el coeficiente es negativo) la probabilidad de que se elija la
alternativa j, en comparación con la categoría base, pero no cuantifican la variación.
Para contar con una medida que cuantifique el impacto de determinada variable sobre la probabilidad de
elegir la alternativa j, se recurre a los efectos marginales, los cuales, para el Modelo Logit Multinomial, están
dados por lo siguiente:
T
e wi γ j
Pji = T T T
ewi γalt1 + ewi γalt2 + ... + ewi γaltm
en donde altj hace referencia a la alternativa j (alt1 se refiere a la alternativa 1, alt2 se refiere a la alternativa
2 y así sucesivamente). Ahora, definiendo wiT γalth como zhi , se tiene:
ez ji
Pji =
ez1i + ez2i + ... + ezmi
Por lo tanto, como zhi = wiT γalth = γh0 + γh1 w1i + ... + γhr wri [el primer subíndice (en los coeficien-
tes) hace referencia a la alternativa a la que corresponden los coeficientes (m alternativas), el segundo subín-
dice de los coeficientes (primer subíndice para las variables) hace referencia a la variable alternative-invariant
(r variables alternative-invariant) y el segundo subíndice de las variables hace referencia al individuo (entidad)],
entonces se tiene que:
γ j f ez ji (∑m zhi z ji
∑m zhi
∂Pji h =1 e ) − e h =1 γh f e
=
(∑m zhi 2
h =1 e )
∂w f
γ j f ez ji (∑m zhi ez ji ∑m zhi
∂Pji h =1 e ) h =1 γh f e
= −
(∑m zhi 2 (∑m zhi 2
h =1 e ) h =1 e )
∂w f
ez ji ∑m h =1 γh f e
zhi
∂Pji
= m z γj f −
∂w f ∑h=1 e hi ∑m h =1 e
zhi
!
∂Pji m
= Pji γ j f − ∑ Phi γh f
∂w f h =1
Debe notarse que los efectos marginales no necesariamente poseen el mismo signo que los coeficientes; por
lo tanto, los coeficientes indican (no miden) si la probabilidad de elegir la alternativa j aumenta o disminuye en
comparación con la categoría base, y los efectos marginales cuantifican la variación en la probabilidad de elegir
determinada alternativa ante la variación en una variable explicativa (alternative-invariant) específica.
Al igual que con los modelos Logit y Probit estudiados en el capítulo anterior, en el Modelo Logit Multino-
mial, para una misma variable se obtienen n efectos marginales distintos (hay n individuos (entidades)). Por lo
tanto, para obtener una medida general del efecto marginal de una variable particular, pueden emplearse los
dos “caminos” expuestos en el Capítulo 6 de esta obra (al lector que no tenga claridad al respecto, se sugiere
consultar dicho capítulo, en cuanto la lógica aplicada al actual es la misma).
Es importante señalar que los efectos marginales de la misma variable deben sumar cero entre las alterna-
tivas, es decir:
102
m ∂Pj
∑ ∂w f =0
j =1
∂P
en donde ∂wjf es el efecto marginal de la variable w f sobre la probabilidad de elegir la alternativa j y exis-
ten m alternativas. Asimismo, es importante traer a colación que sin importar la alternativa que se elija como
categoría base, aunque los coeficientes cambien, los efectos marginales serán los mismos.
Hasta este punto se ha dado una muy breve exposición de los fundamentos teóricos sobre los que se es-
tructura el modelo que ahora pretendemos abordar de manera práctica en Python. Es posible que la estructura
matemática sobre la que se construye el Modelo Logit Multinomial no sea totalmente comprensible para el
lector; sin embargo, puede que la labor práctica le sea de ayuda para tener mayor claridad sobre determinados
aspectos.
El anterior bloque de código permite cargar las herramientas necesarias para llevar a cabo las acciones
requeridas en el análisis econométrico planteado; en caso de que el lector no tenga claridad al respecto, se le
pide consultar la sección “Preparación del entorno” del Capítulo 2, en donde se brinda una explicación sobre
tal cuestión.
103
Una vez ejecutada la línea de código de esta celda, el dataset debería haber sido importado correctamente.
Para comprobarlo, recurrimos al método .head(), aplicado al objeto que contiene el dataset (en este caso le
hemos asignado el nombre data):
In [3]: [Link]()
Podemos observar que el proceso de importación ha sido exitoso y la información del dataset ha sido alma-
cenada correctamente en un DataFrame de pandas. Ahora, podemos continuar tranquilamente con el desarrollo
de nuestro análisis, en cuanto hemos concluido satisfactoriamente el primer paso. Procederemos a construir y
estimar un modelo logit multinomial, el cual toma como variable dependiente una variable discreta no binaria;
por lo tanto, resulta de utilidad conocer las alternativas que contiene dicha variable.
Para conocer las alternativas de la variable dependiente basta con seleccionarla (empleando los corchetes
’[ ]’) y aplicarle el método .unique(), el cual permite conocer los valores únicos que tiene un objeto. El código
a ejecutar es:
In [4]: data["mode"].unique()
El resultado que obtenemos es un ndarray de NumPy, en el cual se nos informa que los valores únicos
(alternativas) que toma la variable dependiente son ’charter’, ’private’, ’pier’ y ’beach’; por lo tanto, contamos
con cuatro alternativas y, consecuentemente, obtendremos tres sets de coeficientes (el lector ya debería saber
el porqué).
Para nuestro caso es muy sencillo saber el número de alternativas con las que se cuenta, pues basta con
enumerarlas; sin embargo, si la variable tuviera muchas alternativas resultaría tedioso contarlas una a una.
pandas ofrece una solución efectiva a dicho inconveniente: en vez de aplicar el método .unique() se utiliza el
método .nunique(), el cual cuenta el número de valores únicos. Para comprobar que efectivamente se tienen
4 alternativas, se ejecuta el siguiente código:
104
In [5]: data["mode"].nunique()
Out[5]: 4
In [6]: Y = data["mode"]
X = data["income"]
Construcción y estimación
Para la construcción y estimación del modelo se emplean la función MNLogit(...) de Statsmodels y el
método .fit(); para la visualización de los resultados, la función print(...) y el método .summary(). El
código a ejecutar es el siguiente:
105
mode=private coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
const 0.7389 0.197 3.756 0.000 0.353 1.125
income 0.0919 0.041 2.260 0.024 0.012 0.172
================================================================================
Nota: La información con la que trabaja la función MNLogit(...) de Statsmodels debe estar contenida
en un dataset de formato wide. En caso de contar con un formato long, el dataset debe reorganizarse de
modo que, antes de emplear su información en la función MNLogit(...) de Statsmodels, su formato sea
wide.
El lector debe recordar que se mencionó que, para la estimación, el set de coeficientes de una alternativa
se normaliza a cero y que la interpretación de los coeficientes obtenidos se realiza con referencia a la categoría
base, sin embargo, ¿cuál es la categoría base?
Para este ejemplo concreto, la categoría base es la alternativa ’beach’, pero ¿por qué? La respuesta se en-
cuentra en la lógica que aplica la función MNLogit(...) de Statsmodels: cuando se pasa como variable depen-
diente una variable no codificada por números (el lector debe observar que la variable mode no está codificada),
la función ejecuta la codificación automáticamente de forma alfabética y toma como categoría base la alterna-
tiva que ocupa la primera posición en dicho orden. Así, como las alternativas de la variable mode son ’charter’,
’private’, ’beach’ y ’pier’, la categoría base es ’beach’, pues es la que ocupa la primera posición si se organizan
alfabéticamente tales alternativas.
Ahora, la interpretación de los coeficientes se debe realizar con respecto a la alternativa ’beach’. Así, se tiene
que en comparación con pesca en la playa (’beach’), un ingreso más alto está asociado a una menor probabilidad
de pesca en bote alquilado (’charter’) o en el muelle (’pier’) y a una mayor probabilidad de pesca en bote privado
(’private’).
Hemos dicho que un mayor ingreso está asociado, en comparación con pesca en la playa, a una probabilidad
mayor o menor de determinadas alternativas, pero ¿de cuánto es esta mayor o menor probabilidad? Recuerde
que los coeficientes no cuantifican la relación, tan solo indican su sentido; para obtener una medida del im-
pacto de determinada variable sobre la probabilidad de elegir una alternativa específica se emplean los efectos
marginales.
Para obtener los efectos marginales de un modelo logit multinomial basta con aplicar el método
.get_margeff() a la instancia de resultados de dicho modelo; para visualizar estos efectos marginales, tan
solo se requiere el empleo del método .summary() y de la función print(...), así:
106
--------------------------------------------------------------------------------
income -0.0112 0.006 -1.876 0.061 -0.023 0.000
--------------------------------------------------------------------------------
mode=pier dy/dx std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
income -0.0208 0.005 -4.040 0.000 -0.031 -0.011
------------------------------------------------------------------------------
mode=private dy/dx std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
income 0.0318 0.005 6.039 0.000 0.021 0.042
================================================================================
De este modo, se tiene que un incremento de una unidad en el ingreso está asociado a que la pesca en la
playa (’beach’) sea 0.02 puntos porcentuales más probable; a que la pesca en bote alquilado (’charter’) sea 1.12
puntos porcentuales menos probable; a que la pesca en el muelle (’pier’) sea 2.08 puntos porcentuales menos
probable; y a que la pesca en bote privado (’private’) sea 3.18 puntos porcentuales más probable. (El lector
puede comprobar que la suma de estos efectos marginales es cero)
El lector debe tener en cuenta que existen varias maneras de obtener una medida general del efec-
to marginal de determinada variable; los efectos marginales que son calculados por defecto por el método
.get_margeff() son los efectos marginales promedio. Si lo que se quiere conocer son los efectos marginales
en la media, se debe asignar el valor de ’mean’ al parámetro at del método .get_margeff(), así:
107
Ahora, se tiene que un incremento de una unidad en el ingreso está asociado a que la pesca en la playa
(’beach’) sea 0.008 p.p más probable; a que la pesca en bote alquilado (’charter’) sea 1.2 p.p menos probable;
a que la pesca en el muelle (’pier’) sea 2.07 p.p menos probable; y a que la pesca en bote privado (’private’) sea
3.26 p.p más probable. (Nuevamente, el lector puede comprobar que la suma de estos efectos marginales es
cero).
Si el lector contrasta los resultados obtenidos con los de Katchova (2013, 3) notará que los que ella obtiene
son diferentes; esto se debe a que en el ejercicio desarrollado por la profesora Katchova no se toma como
categoría base ’beach’ sino ’charter’, ¿es posible que nosotros hagamos lo mismo? La respuesta es sí.
La función MNLogit(...) de Statsmodels toma como categoría base la alternativa que ocupa la primera
posición en orden alfabético o numérico y no es posible modificar tal lógica en su funcionamiento. Sin embargo,
si se recodifican las alternativas de modo que la que se desee tener como categoría base ocupe la primera
posición, entonces la función MNLogit(...) de Statsmodels sí asumirá dicha alternativa como la categoría
base.
A modo de ejemplo, recurriendo a herramientas de pandas, se puede codificar una variable categórica usando
el método .map(...) y empleando un diccionario Python “{...}” como argumento; dentro de dicho diccionario
las claves corresponderán a las alternativas y los valores a los códigos asignados a cada una.
Por ejemplo, a continuación se crea una variable ’mode2’ que corresponde a la variable ’mode’ codificada
de modo que la alternativa ’charter’ se encuentre en la primera posición en orden numérico y, por ende, sea
tomada como categoría base por la función MNLogit(...) de Statsmodels:
Empleando la misma estructura de código que hemos utilizado con anterioridad, se obtiene el siguiente
modelo:
In [11]: Y = data["mode2"]
ModeloLogitMN = [Link](Y, sm.add_constant(X))
ResultadosLogitMN = [Link]()
print([Link]())
108
income 0.0316 0.042 0.756 0.450 -0.050 0.114
-----------------------------------------------------------------------------------
mode2=3_pier coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------
const -0.5271 0.178 -2.965 0.003 -0.876 -0.179
income -0.1118 0.044 -2.541 0.011 -0.198 -0.026
--------------------------------------------------------------------------------
mode2=4_private coef std err z P>|z| [0.025 0.975]
-----------------------------------------------------------------------------------
const -0.6024 0.136 -4.426 0.000 -0.869 -0.336
income 0.1235 0.028 4.426 0.000 0.069 0.178
===================================================================================
Como se puede apreciar, las magnitudes de los coeficientes han cambiado y, además, ahora su interpretación
no se debe hacer con referencia a ’beach’ sino a ’charter’. Así, se tiene que, en comparación con pesca en bote
alquilado (’charter’), un ingreso más alto está asociado a una mayor probabilidad de pesca en la playa (’beach’)
y en bote privado (’private’) y a una menor probabilidad de pesca en el muelle (’pier’). ¿Qué hay de los efectos
marginales?
El lector debe recordar que, sin importar cuál sea la categoría base, el efecto marginal de una variable
específica, a diferencia del coeficiente, siempre será el mismo. Para comprobarlo, basta con aplicar el método
.get_margeff() a la instancia de resultados del modelo:
109
Se puede observar que los efectos marginales obtenidos del modelo con categoría base ’charter’ son exac-
tamente iguales a los obtenidos del modelo con categoría base ’beach’ (puede examinarlos uno por uno si así
lo desea). Recordando que el método .get_margeff() calcula por defecto efectos marginales promedio, si lo
que se desea conocer son los efectos marginales en la media, tan solo se requiere asignar el valor de ’mean’ al
parámetro at:
Nuevamente, los efectos marginales en la media son exactamente iguales para el modelo con categoría base
’beach’ y para el modelo con categoría base ’charter’, igualdad que se mantiene para cualquier otra alternativa
de la variable dependiente (empleando, obviamente, los mismos datos).
En este punto se cierra este capítulo, con el cual se concluye esta breve obra, la cual ha buscado ofrecer
ejemplos ilustrativos de análisis econométrico de nivel básico con el lenguaje de programación Python y en la
que se han presentado algunas de las herramientas más relevantes para dicho ejercicio, tratando de brindar al
lector breves, pero útiles, descripciones del funcionamiento de cada una de estas.
110
Referencias
Chatterjee, S., y Simonoff, J. (2013). Handbook of Regression Analysis. Estados Unidos: John Wiley & Sons,
Inc.
Dormann, C.F., Elith, J., Bacher, S., Buchmann, C., Gudrun, C., Carré, G., García, J.R., Gruber, B., Lafourcade,
B., Leitão, P.J., Münkemüller, T., McClean, C., Osborne, P.E., Reineking, B., Schröder, B., Skidmore, A.K.,
Zurell, D., y Lautenbach, S. (2013). Collinearity: a review of methods to deal with it and a simulation
study evaluating their performance. Ecograpghy, 36, 27-46. doi: 10.1111/j.1600-0587.2012.07348.x
Fahrmeier, L., Kneib, T., y Lang, Stefan. (2007). Regression. Modelle, Methoden und Anwendungen. Springer.
Greene, W.H. (2003). Econometric Analysis. New Jersey, Estados Unidos: Pearson Education, Inc.
Gujarati, D.N. y Porter, D.C. (2010). Econometría. México: McGraw-Hill/Interamericana Editores, S.A. de
C.V.
James, G., Witten, D., Hastie, T., y Tibshirani, R. (2013). An Introduction to Statistical Learning with Applica-
tions in R. doi: 10.1007/978-1-4614-7138-7
Katchova, A. (2013). Econometrics – Probit and Logit Models. Recuperado de: [Link]
le/d/0BwogTI8d6EEidjg2OGl4b3dxVmc/edit
Katchova, A. (2013,2). Econometrics – Multinomial Probit and Logit Models. Recuperado de:
[Link]
Katchova, A. (2013, 3). Econometrics – Multinomial Probit and Logit Models, Conditional Lo-
git Model, Mixed Logit Model Examples. Recuperado de: [Link]
TI8d6EEiLVZ0N1h3LTBLYlk/edit
Kleiber, C., y Zeileis, A. (2008). Applied Econometrics with R. Estados Unidos: Springer Science+Business
Media, LLC.
Murray, L., Nguyen, H., Lee, Y-F., Remmenga, M., y Smith, D.W. (2012). Variance Inflation Factors in Re-
gression Models with dummy variables. Conference on Applied Statistics in Agriculture. doi: 10.4148/2475-
7772.1034
Perktold, J., Seabold, Skipper., y Taylor, J. (2018). Welcome to Statsmodels’s Documentation. Recuperado
de: [Link]
111
Python Software Foundation. (s.f.). What is Python? Executive Summary. Recuperado de: [Link]
[Link]/doc/essays/blurb/
Vaart, A.W. van der. (1998). Asymptotic Statistics. Reino Unido: Cambridge University Press.
Verbeek, M. (2004). A guide to modern econometrics. Inglaterra: John Wiley & Sons Ltd.
Wooldridge, J.M. (2010). Introducción a la econometría. Un enfoque moderno. México: Cengage Learning
Editores, S.A. de C.V.
Würtz, D., y Katzgraber, H.G. (2009). Precise finite-sample quantiles of the Jarque-Bera adjusted Lagran-
ge multiplier test. ETH Econohysics Working and White Papers Series. Recuperado de: [Link]
[Link]/sites/default/files/[Link]
112