Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Unprecise p-values in Stargazer

Tags:

r

stargazer

I want the same stars for significancies in regression output in stargazer as in the "normal output".

I produce data

library("stargazer"); library("lmtest"); library("sandwich")
set.seed(1234)
df <- data.frame(y=1001:1100)
df$x <- c(1:70,-100:-71) + rnorm(100, 0, 74.8)
model <- lm(log(y) ~ x, data=df)

and get some model estimates where the coefficient on x has a p-value of 0.1023

coeftest(model, vcov = vcovHC(model, type="HC3"))

I want to have these results in LaTeX. Based on the same function I calculate heteroscedasticity consistent standard estimates and let stargazer use them.

stderr_HC3_model <- sqrt(diag(vcovHC(model, type = "HC3")))
stargazer(model, se=list(stderr_HC3_model))

The stargazer output has a star at the coefficient indicating significance when alpha=10%. I want stargazer to give the same as the coeftest. (Because of the comparability with Stata where reg L_y x, vce(hc3) gives exactly the coeftest results.)

I played around with the stargazer options p.auto, t.auto which did not help. When I execute "stargazer" I cannot view the underlying code as it is possible in other cases. What to do?


Richards answer helped me. I indicate the steps I used to give out more than one regression (let's say ols_a and ols_b).

ses <- list(coeftest(ols_a, vcov = vcovHC(ols_a, type="HC3"))[,2],
        coeftest(ols_b, vcov = vcovHC(ols_b, type="HC3"))[,2])
pvals <- list(coeftest(ols_a, vcov = vcovHC(ols_a, type="HC3"))[,4],
          coeftest(ols_b, vcov = vcovHC(ols_b, type="HC3"))[,4])
stargazer(ols_a, ols_b, type="text", p=pvals, se=ses)
like image 552
user3070843 Avatar asked Apr 22 '14 11:04

user3070843


3 Answers

There are inherent dangers associated with the se argument.

When using this approach, the user should be cautious wrt the arguments t.auto and p.auto, both of which default to TRUE. I think it would be cautious to set them both to FALSE, and supply manually t and p values.

Failure to do so, and you risk getting significance stars not in sync with the displayed p-values. (I suspect that stargazer will simply reuse the se, which are now different from the default ones, and recompute the displayed stars using this input; which will naturally yield unexpected results.)

See also:

  • Displaying p-values instead of SEs in parenthesis
like image 151
landroni Avatar answered Nov 02 '22 10:11

landroni


You need to provide the p values associated with your coeftest. From the man page.

p a list of numeric vectors that will replace the default p-values for each model. Matched by element names. These will form the basis of decisions about significance stars

The following should work.

test <- coeftest(model, vcov = vcovHC(model, type="HC3"))
ses <- test[, 2]
pvals <- test[, 4]
stargazer(model, type="text", p=pvals, se=ses)

This provides the following.

===============================================
                        Dependent variable:    
                    ---------------------------
                              log(y)           
-----------------------------------------------
x                            -0.00005          


Constant                     6.956***          
                              (0.003)          

-----------------------------------------------
Observations                    100            
R2                             0.026           
Adjusted R2                    0.016           
Residual Std. Error       0.027 (df = 98)      
F Statistic             2.620 (df = 1; 98)     
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01
like image 32
Richard Herron Avatar answered Nov 02 '22 10:11

Richard Herron


It may be a minor issue but Richard's answer is actually not entirely correct - his stargazer output does not report any standard errors nor potential significance stars for the variable x.

Also when reporting only a single model in stargazer manual coefficients, se, p and t values have to be provided in a list. Otherwise stargazer will report an empty list.

The (slightly) corrected example:

test <- coeftest(model, vcov = vcovHC(model, type="HC3"))
ses <- list(test[, 2])
pvals <- list(test[, 4])
stargazer(model, type="text", p=pvals, se=ses)

Output:

=======================================================================
                                         Dependent variable:           
                              -----------------------------------------
                                        Daily added investors          
                                              negative                 
                                              binomial                 
-----------------------------------------------------------------------
log(lag_raised_amount + 1)                    -0.466***                
                                               (0.124)                 

lag_target1                                   -0.661***                
                                               (0.134)                 

Constant                                      -3.480**                 
                                               (1.290)                 

-----------------------------------------------------------------------
Observations                                    6,513                  
Log Likelihood                                 -8,834                
theta                                     1.840*** (0.081)             
Akaike Inf. Crit.                              17,924                
=======================================================================
Note:                         + p<0.1; * p<0.05; ** p<0.01; *** p<0.001
like image 26
JNWHH Avatar answered Nov 02 '22 09:11

JNWHH