Un ejemplo sobre cómo combinar Stata y R para análisis de datos

Ya en algún “post” anterior había dicho que eso de hacer análisis de datos en R no me gusta mucho. Poco a poco le he ido agarrando más el gusto; sin embargo, Stata sigue siendo mi herramienta favorita para análisis de datos (por cierto, para los “fans” de Stata, hace unos días anunciaron su versión número 13 -si tienen algún dinerito extra ya saben qué pueden hacer con él). ¿Por qué no combinar ambas herramientas? ¿Qué tal si usamos Stata para análisis de datos y R para hacer gráficas bonitas? La idea suena bien, especialmente si no hay que codificar mucho para obtener los resultados que queremos.

Aquí va un ejemplo. Si quieren seguir los mismos pasos del ejemplo pueden descargar la base de datos que voy a utilizar dando click aquí (es una base de datos con números aleatorios que servirá para los propósitos de este post). Una vez descargada la base, pongan el archivo en la carpeta que deseen y utilicen esa misma carpeta como el “working directory” de la sesión de Stata y R que vayan a comenzar. Mi carpeta se llama “StataR” y se encuentra en mi “Desktop”. Así, para Stata usé cd "~/Desktop/StataR" y para R usé setwd("~/Desktop/StataR").

La base del ejemplo es un panel con siete paneles (valga la redundancia) y 10 observaciones por panel. El objetivo será usar los siguientes estimadores y comparar sus estimados gráficamente: 1) una regresión simple de mínimos cuadrados ordinarios, 2) una regresión de efectos fijos, y 3) una regresión de efectos aleatorios. Antes de empezar hay que instalar en Stata el módulo “svmatf” para convertir matrices en archivos con terminación .dta (aquí la información sobre svmatf). También, antes que otra cosa hay que “declarar” que la base es de tipo panel: xtset country year. (NOTA:  escribí este post con la versión 12 de Stata, la cual aún no tenía el comando putexcel. Este comando podría sustituir a svmatf para crear la o las bases de datos requeridas. Por acá les dejo más información sobre putexcel).

Empecemos, pues. Lo primero es crear un “loop” que haga todo automáticamente. Línea por línea el loop indica lo siguiente:

  1. foreach var of varlist y {  inicia el loop indicando a Stata que queremos usar la variable “y” en ese loop.
  2. regress `var’ x1 x2 x3   indica a Stata que queremos correr una regresión de mínimos cuadrados ordinarios con nuestra variable y como variable dependiente (ahora denominada `var’ gracias al loop) y x1, x2 y x3 como variables independientes.
  3. estimates store ols_`var’  indica a Stata que queremos guardar los estimados de la regresión anterior.
  4. matrix mat_ols_`var’ = ([_b[x1],_se[x1]]\[_b[x2],_se[x2]]\[_b[x3],_se[x3]])  crea una matriz de tres renglones y dos columnas. La primera columna con los coeficientes y la segunda columna con los errores estándar.
  5. svmatf, mat(mat_ols_`var’) fil(mat_ols_`var’.dta)  indica que queremos guardar la matriz del punto 4 como una base de datos con terminación .dta.
  6. Se repite lo mismo pero, en vez de regress, escribimos xtreg para usar efectos fijos y efectos aleatorios.
  7. Al final no olviden cerrar el loop con un corchete: }

Aquí va pues el programa entero:

foreach var of varlist y {
regress `var' x1 x2 x3
estimates store ols_`var'
matrix mat_ols_`var' = ([_b[x1],_se[x1]]\[_b[x2],_se[x2]]\[_b[x3],_se[x3]])
svmatf, mat(mat_ols_`var') fil(mat_ols_`var'.dta)

xtreg `var' x1 x2 x3, fe
estimates store fe_`var'
matrix mat_fe_`var' = ([_b[x1],_se[x1]]\[_b[x2],_se[x2]]\[_b[x3],_se[x3]])
svmatf, mat(mat_fe_`var') fil(mat_fe_`var'.dta)

xtreg `var' x1 x2 x3, re
estimates store re_`var'
matrix mat_re_`var' = ([_b[x1],_se[x1]]\[_b[x2],_se[x2]]\[_b[x3],_se[x3]])
svmatf, mat(mat_re_`var') fil(mat_re_`var'.dta)
}

Después de correr el programa anterior van a encontrar en su “working directory” tres bases de datos llamadas “mat_ols_y.dta”, “mat_fe_y.dta” y “mat_re_y.dta”. Estas bases de datos contienen, cada una, tres renglones y dos columnas: cada renglón corresponde a cada una de las variables independientes de nuestros modelos (x1, x2 y x3 en este caso), la primer columna con información de los coeficientes y la segunda con los errores estándar. Pasemos ahora a R.

Lo primero que queremos hacer es juntar las tres bases anteriores en una sola base de datos. En vez de decirle a R tres veces que haga las cosas, le podemos dar instrucciones una sola vez con el siguiente comando:

archivos = list.files("~/Desktop/StataR/")  #Lee el nombre de los archivos en la carpeta señalada.
archivos = archivos[-4]  #Quita el elemento cuatro que no queremos usar (este elemento es la base panel que utilizamos en Stata: "PanelPrueba.dta").
archivos = paste("~/Desktop/StataR/",archivos,sep="") #Pega el nombre de los archivos anteriores a la carpeta donde tengo estos archivos (así R reconoce dónde se encuentran. Si tienen esta carpeta como working directory esto no es necesario)
datos = do.call("rbind", lapply(archivos, read.dta)) #Junta las tres bases de datos (mat_ols_y.dta, mat_fe_y.dta, etc.) en una sola matriz. Ahora vamos a tener 9 renglones y dos columnas.

El resultado que se obtiene de lo anterior debe ser algo así:
> datos
c1       c2        row
1 2424529152 1156516224 r1
2 1822699648 2028055936 r2
3 309718336 368552480   r3
11 559063936 915594368 r1
21 87445760 349461184   r2
31 92622344 293654976   r3
12 1734069376 999303744 r1
22 420030528 497228768   r2
32 226606960 324458560   r3

Ahora, tengan cuidado pues R leyó los archivos en orden alfabético: primero “mat_fe_y.dta”, posteriormente “mat_ols_y.dta” y al último “mat_re_y.dta”. Así, los tres primeros renglones corresponden a los estimados del modelo de efectos fijos, los renglones 4 a 6 al modelo de mínimos cuadrados ordinarios y los últimos al modelo de efectos aleatorios.
Para identificar bien a que modelo pertenecen, creemos un factor y añadámoslo como cuarta columna a nuestra base “datos”:

datos$modelo = factor(c(rep("FE",3),rep("OLS",3),rep("RE",3)))  #Aquí sólo le estoy diciendo a R que añada a mi base datos un factor de nueve elementos, los primeros tres llamados FE, los siguientes tres OLS y los últimos tres RE.

Ahora sí, estamos listos para hacer la tan ansiada gráfica para comparar visualmente los estimados de los tres modelos. Aquí va el programa con notas aclaratorias en cada renglón (necesitan cargar el paquete “ggplot2″):

rango = qnorm(1 - 0.1 / 2) #Para crear los intervalos de confianza al 90% (si quieren al 95% sustituyan el 0.1 por un 0.05)
grafica = qplot(x=row, y=c1, facets= modelo~., data = datos,  #Indicamos que queremos el nombre de las variables en el eje de las X y el de los coeficientes (llamados aquí c1) en el eje de las Y. Utilizamos un facet para agrupar según el tipo de modelo (usamos el factor que creamos anteriormente llamado modelo)
 ymin = c1 - rango * c2, #Creamos las lineas de los intervalos de confianza (recuerden que c1 son los coeficientes y c2 sus errores estándar.
 ymax = c1 + rango * c2,
 geom = "pointrange",  #Indicamos a ggplot2 que lo que tenemos son estimados con rangos, etc.
 xlab="estimado", #Indica el título del eje de las X.
 ylab="", #Ningún título para este eje pues sólo son los valores que toman los estimados.
 main = "Comparativo modelos \n (IC: 90%) \n")  #Título de la gráfica.

grafica = grafica + coord_flip() #Giramos la gráfica para que ahora el eje de las X pase al eje de las Y (así la comparación será visualmente más sencilla).
grafica = grafica + theme_bw() + theme(axis.text.y = element_text(size=14),
        axis.text.x = element_text(size=13))  #Todo esto es para cuestiones estéticas, tamaño de letra, etc.
grafica = grafica + geom_hline(yintercept = 0, colour = I("#FF0000"), alpha = I(1/2), size=1.2) #ponemos una línea horizontal en 0.
grafica = grafica + geom_point(aes(size=2)) #tamaño de los puntos.
grafica = grafica + theme(legend.position="none") #Quitamos la leyenda que aparece a la derecha.
print(grafica)

Si todo sale bien deben obtener una gráfica igual a la de aquí abajo (r1, r2 y r3 son los nombres que Stata dio a nuestras variables independientes cuando creamos las matrices. Así que r1 corresponde a la variable x1, r2 a la variable x2, etc.). Yupiiii… ojalá este post les sea útil.

ComparativoModelos

About these ads

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: