mtable el equivalente en R de estimates table de stata

Indicaciones de stata
Antes de convertirme en asiduo usuario de R  fui usuario de STATA que satisface bien las necesidades de cálculo siempre y cuando no se tenga que programar, una vez que se requiere programar se tienen dos alternativas:

1. Acostumbrarse a su incómoda sintaxis ('lenguaje' para ser un poco más formal) o
2. Emigrar hacia otro programa cuyo lenguaje sea más manejable, en mi caso ese lenguaje fue R, pero este post no trata sobre las razones por las cuales cambié de programa, sino que se trata de mostrar algo que en stata me gustaba mucho y que pensaba que lo echaría en falta al cambiarme a R, se trata de la comparación de dos (o más) modelos de regresión cuyo output se vea en una sola tabla, en stata esto se logra con la función estimates table y en R con la función mtable del paquete memisc.


Ejemplos de ambos resultados se muestran a continuación

Indicaciones de stata


Los siguientes comandos de stata son un extracto de un do file
  sysuse auto                                  /* Cargar la base de datos                       */
  regress price weight length mpg              /* Estimación una regresión simple con constante */
  estimates store constant                     /* Guardando los resultados de la estimación     */
  regress price weight length mpg , noconstant /* Estimación una regresión simple con constante */
  estimates store nonconstant                  /* Guardando los resultados de la estimación     */

  /* A continuación se crea la tabla para comparar los dos modelos */ 
  estimates table constant nonconstant, b(%9.4f) star stats( r2, r2_a, F, ll, N)
En seguida se ve que el resultado es el siguiente:
--------------------------------------------
    Variable |   constant     nonconstant   
-------------+------------------------------
      weight |    4.3648***      3.2291**   
      length | -104.8682*      -25.4414     
         mpg |  -86.7893        54.6221     
       _cons |  1.45e+04*                   
-------------+------------------------------
          r2 |    0.3574         0.8713     
        r2_a |    0.3298         0.8659     
           F |   12.9762       160.2632     
          ll | -679.3516      -682.4405     
           N |        74             74     
--------------------------------------------
    legend: * p<0.05; ** p<0.01; *** p<0.001

Los argumentos de la función estimates table indican lo siguiente: b(%9.4f)es para redondear al cuarto decimal, star es para indicar que la siginificancia se reflejará con el uso de estrellas (ver leyenda de la tabla), stats es la lista de estadísticos que queremos que aparezcan en la tabla, además de los coeficientes, los estadíscos que se han pedido fueron: R cuadrado (r2), R cuadrado ajustado (r2_a), el estadístico F (F), log-likelihood (ll) y tamaño muestral (N), se pueden agregar otros más para igualar el resultado de la tabla de R.

El resultado es muy elegante y práctico, sin embargo tiene unas limitantes, sólo se indica el nivel de significancia a partir de las 'estrellas' pero no se pueden proporcionar los valores de los errores estándar, eso lo dice claro en el manual de ayuda (help(estimates table)) cuando dice:“The star and star() options may not be combined with the se, t, or p options”. Lo cual es una verdadera lástima.

Un resultado más limpio se podría lograr con la función .esttab del paquete estout, sin embargo, este no lo ilustro porque es más común usar estimates store.

Indicaciones de R


library(foreign)
auto <- read.dta("http://www.stata-press.com/data/r9/auto.dta")  # Importar datos
const <- lm(price ~ weight + length + mpg, data = auto)  # Modelo con constante
nonconst <- lm(price ~ weight + length + mpg - 1, data = auto)  # Modelo sin constante
library(memisc)
mtable(const, nonconst)
## 
## Calls:
## const: lm(formula = price ~ weight + length + mpg, data = auto)
## nonconst: lm(formula = price ~ weight + length + mpg - 1, data = auto)
## 
## ===========================================
##                     const       nonconst   
## -------------------------------------------
## (Intercept)     14542.434*                 
##                 (5890.632)                 
## weight              4.365***     3.229**   
##                    (1.167)      (1.111)    
## length           -104.868*     -25.441     
##                   (39.722)     (24.117)    
## mpg               -86.789       54.622     
##                   (83.943)     (63.526)    
## -------------------------------------------
## R-squared               0.357         0.871
## adj. R-squared          0.330         0.866
## sigma                2414.563      2499.692
## F                      12.976       160.263
## p                       0.000         0.000
## Log-likelihood       -679.352      -682.440
## Deviance        408107984.249 443640629.333
## AIC                  1368.703      1372.881
## BIC                  1380.224      1382.097
## N                      74            74    
## ===========================================

Con lo cual cambiarse de stata a R no es tan traumático.

Resultados en formato latex también es posible tanto en stata como en R, pero eso es tema para otra entrada.


4 comentarios:

  1. Buenos días.
    Muy interesante tu blog sobre R. Tengo un problema con la libreria "memics". No encuentra dicha librería. Pongo el comando y el mensaje de error:
    > library(memisc)
    Error in library(memisc) : there is no package called ‘memisc’
    Alguien tiene alguna idea de por qué no funciona.
    José Julián Escario
    Huesca
    España

    ResponderEliminar
  2. Eso es porque el paquete 'memisc' no está instalado. Debés instalarlo primero con 'install.packages("memisc")', una vez instalado, lo podés cargar en tu sesión de trabajo vía 'library(memisc)'.

    ResponderEliminar
  3. Hola, buenos días.

    Muy buen explicación. Pero me surge una duda. ¿Cómo haces para que en la tabla salgan todos los criterios de información (AIC BIC Devian) y demás criterios? Si pongo mtable(model1, model2) me aparece solo hasta aquí

    ## ===========================================
    ## model 1 model 2
    ## -------------------------------------------
    ## (Intercept) 14542.434*
    ## (5890.632)
    ## weight 4.365*** 3.229**
    ## (1.167) (1.111)
    ## length -104.868* -25.441
    ## (39.722) (24.117)
    ## mpg -86.789 54.622
    ## (83.943) (63.526)
    ## -------------------------------------------

    Qué código utilizas para agregar esos criterios a la tabla?

    Muchas gracias.

    ResponderEliminar
    Respuestas
    1. Los códigos publicados fueron los utilizados para obtener ese output. Sin embargo, probá con el siguiente código `mtable(model1, model2, summary.stats=TRUE)`

      Eliminar